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

r465

Browse Source
This commit is contained in:
chhylp123
2022-12-02 09:40:25 -05:00
parent a4d6cebc32
commit 786f8bdff0
6 changed files with 967 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.h
+1 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.17.7-r461"
#define HA_VERSION "0.18.0-r465"
#define VERBOSE 0
Overlaps.cpp
+160 -7
View File
@@ -7718,7 +7718,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen)
kv_push(uint32_t, b_r, w);
aw = asg_arc_a(g, w);
min_edge = (u_int32_t)-1;
min_edge = (uint32_t)-1;
for (t = 0; t < nw; t++)
{
if(aw[t].del) continue;
@@ -13463,6 +13463,153 @@ void hic_clean(asg_t* read_g)
kv_destroy(ax);
}
void hic_clean_adv(asg_t *sg, ug_opt_t *uopt)
{
uint32_t i, k, m, z, v, w; ma_utg_t *u = NULL; uint32_t *ba, bn, n_vtx, beg, end, n0, n1;
ma_ug_t *ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); n_vtx = ug->g->n_seq<<1; double bub_rate = 0.1;
uint8_t *bf = NULL; bubble_type *bub = gen_bubble_chain(sg, ug, uopt, &bf);
uint64_t tLen, vocc, socc, pocc; buf_t b; memset(&b, 0, sizeof(buf_t)); CALLOC(b.a, n_vtx);
REALLOC(bf, n_vtx); memset(bf, 0, sizeof((*bf))*n_vtx);
kvec_t(uint64_t) buf; kv_init(buf); n0 = n1 = 0;
for (i = 0; i < bub->b_ug->u.n; i++) {
u = &(bub->b_ug->u.a[i]);
if(u->n == 0) continue;
for (k = 0; k < u->n; k++) {///bubble chain
get_bubbles(bub, u->a[k]>>33, &beg, &end, &ba, &bn, NULL);///bubble
for (m = vocc = tLen = 0; m < bn; m++) {
bf[ba[m]] = bf[ba[m]^1] = 1;
tLen += ug->u.a[ba[m]>>1].len;
if(IF_HOM((ba[m]>>1), *bub)) continue;
if(ug->g->seq[ba[m]>>1].del) continue;
vocc += ug->u.a[ba[m]>>1].n;
}
if(beg != (uint32_t)-1) tLen += ug->u.a[beg>>1].len;
if(end != (uint32_t)-1) tLen += ug->u.a[end>>1].len;
if(vocc) {
for (m = buf.n = 0; m < bn; m++) {
v = ba[m];
if(ug->g->seq[v>>1].del) continue;
if(asg_arc_n(ug->g, v) < 2) continue;
if(get_real_length(ug->g, v, NULL) < 2) continue;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) {
//beg is v, end is b.S.a[0]
//note b.b include end, does not include beg
for (z = socc = 0; z < b.b.n; z++) {
if(b.b.a[z]==v || b.b.a[z]==b.S.a[0]) continue;
socc += ug->u.a[b.b.a[z]>>1].n;
if((!bf[b.b.a[z]])&&(!bf[b.b.a[z]^1])) break;
}
if(z < b.b.n) continue;
kv_push(uint64_t, buf, ((socc<<32)|v));
}
v ^= 1;
if(asg_arc_n(ug->g, v) < 2) continue;
if(get_real_length(ug->g, v, NULL) < 2) continue;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) {
//beg is v, end is b.S.a[0]
//note b.b include end, does not include beg
for (z = socc = 0; z < b.b.n; z++) {
if(b.b.a[z]==v || b.b.a[z]==b.S.a[0]) continue;
socc += ug->u.a[b.b.a[z]>>1].n;
if((!bf[b.b.a[z]])&&(!bf[b.b.a[z]^1])) break;
}
if(z < b.b.n) continue;
kv_push(uint64_t, buf, ((socc<<32)|v));
}
}
radix_sort_arch64(buf.a, buf.a + buf.n);
for (m = pocc = 0; m < buf.n; m++) {
v = (uint32_t)buf.a[m];
if(ug->g->seq[v>>1].del) continue;
if(asg_arc_n(ug->g, v) < 2) continue;
if(get_real_length(ug->g, v, NULL) < 2) continue;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) {
//beg is v, end is b.S.a[0]
//note b.b include end, does not include beg
for (z = socc = 0; z < b.b.n; z++) {
if(b.b.a[z]==v || b.b.a[z]==b.S.a[0]) continue;
socc += ug->u.a[b.b.a[z]>>1].n;
if((!bf[b.b.a[z]])&&(!bf[b.b.a[z]^1])) break;
}
if(z < b.b.n) continue;
if((pocc+socc) >= (vocc*bub_rate)) continue;
pocc += socc;
asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0, NULL);
// fprintf(stderr, "+utg%.6dl\tutg%.6dl\n", (int32_t)(v>>1)+1, (int32_t)(b.S.a[0]>>1)+1);
n0++;
}
}
}
for (m = 0; m < bn; m++) {
bf[ba[m]] = bf[ba[m]^1] = 0;
}
}
}
tLen = get_bub_pop_max_dist_advance(ug->g, &b);
for (v = buf.n = 0; v < n_vtx; ++v) {
if(ug->g->seq[v>>1].del) continue;
if(asg_arc_n(ug->g, v) < 2) continue;
if(get_real_length(ug->g, v, NULL) < 2) continue;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) {
//beg is v, end is b.S.a[0]
//note b.b include end, does not include beg
for (i = socc = 0; i < b.b.n; i++) {
if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue;
socc += ug->u.a[b.b.a[i]>>1].n;
}
if(socc <= 16) kv_push(uint64_t, buf, ((socc<<32)|v));
}
}
uint32_t convex; long long ll, tmp, max_stop_nodeLen, max_stop_baseLen;
radix_sort_arch64(buf.a, buf.a + buf.n);
for (m = 0; m < buf.n; m++) {
v = (uint32_t)buf.a[m];
if(ug->g->seq[v>>1].del) continue;
if(asg_arc_n(ug->g, v) < 2) continue;
if(get_real_length(ug->g, v, NULL) < 2) continue;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) {
w = b.S.a[0]^1;
for (i = socc = 0; i < b.b.n; i++) {
if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue;
socc += ug->u.a[b.b.a[i]>>1].n;
}
if(socc <= 16) {
b.b.n = 0;
get_unitig(ug->g, NULL, v^1, &convex, &ll, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, &b);
for (k = vocc = 0; k < b.b.n; k++) {
if(IF_HOM((b.b.a[k]>>1), *bub)) break;
vocc += ug->u.a[b.b.a[k]>>1].n;
}
if((socc) >= (vocc*bub_rate)) continue;
b.b.n = 0;
get_unitig(ug->g, NULL, w^1, &convex, &ll, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, &b);
for (k = vocc = 0; k < b.b.n; k++) {
if(IF_HOM((b.b.a[k]>>1), *bub)) break;
vocc += ug->u.a[b.b.a[k]>>1].n;
}
if((socc) >= (vocc*bub_rate)) continue;
asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0, NULL);
// fprintf(stderr, "-utg%.6dl\tutg%.6dl\n", (int32_t)(v>>1)+1, (int32_t)(b.S.a[0]>>1)+1);
n1++;
}
}
}
filter_sg_by_ug(sg, ug, uopt);
free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
ma_ug_destroy(ug); free(bf); kv_destroy(buf);
destory_bubbles(bub); free(bub);
// fprintf(stderr, "[M::%s::] # type0::%u, # type1::%u\n", __func__, n0, n1);
}
void update_dump_trio(uint8_t* trio_flag, uint32_t rn, uint8_t *rf, ma_ug_t *ug)
{
uint32_t k, i, x;
@@ -18387,7 +18534,7 @@ void chain_origin_trans_uid_s_bubble(buf_t *pri, buf_t* aux, uint32_t beg, uint3
if(av[i].v == pri_v) priEnd = ((pri_len > av[i].ol)? (pri_len - av[i].ol - 1) : 0);
if(av[i].v == aux_v) auxEnd = ((aux_len > av[i].ol)? (aux_len - av[i].ol - 1) : 0);
}
///[priBeg, priEnd) && [auxBeg, auxEnd)
if(priBeg == (uint32_t)-1 || priEnd == (uint32_t)-1 || auxBeg == (uint32_t)-1 || auxEnd == (uint32_t)-1)
{
fprintf(stderr, "ERROR-s_bubble\n");
@@ -18508,7 +18655,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
for (k = 0; k < p->n; k++)
{
rId = p->a[k]>>33;
cov->cov[rId] += (uCov * cov->read_g->seq[rId].len);
cov->cov[rId] += (uCov * cov->read_g->seq[rId].len);///this the average coverage of the whole bubble
}
}
v = u;
@@ -18517,6 +18664,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
if(t_ch)
{
///is this requirement appropriate?
if(get_real_length(ug->g, v0, NULL) == 2 && get_real_length(ug->g, b->S.a[0]^1, NULL) == 2)
{
long long tmp, max_stop_nodeLen, max_stop_baseLen, bch_occ[2];
@@ -18530,16 +18678,17 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
&max_stop_nodeLen, &max_stop_baseLen, 1, NULL);
get_unitig(ug->g, NULL, bch[1], &convex[1], &bch_occ[1], &tmp,
&max_stop_nodeLen, &max_stop_baseLen, 1, NULL);
///if this is a simple bubble
if(((bch_occ[0] + bch_occ[1] + 1) == (uint32_t)b->b.n) &&
get_real_length(ug->g, convex[0], NULL) == 1 && get_real_length(ug->g, convex[1], NULL) == 1)
{
get_real_length(ug->g, convex[0], &convex[0]);
get_real_length(ug->g, convex[1], &convex[1]);
if(convex[0] == b->S.a[0] && convex[1] == b->S.a[0])
if(convex[0] == b->S.a[0] && convex[1] == b->S.a[0])///double check if it is a simple bubble
{
t_ch->b_buf_0.b.n = 0;
get_unitig(ug->g, NULL, bch[0], &convex[0], &bch_occ[0], &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, &(t_ch->b_buf_0));
for (i = 0; i < t_ch->b_buf_0.b.n; ++i)
for (i = 0; i < t_ch->b_buf_0.b.n; ++i)///retrive one side of the bubble
{
uId = t_ch->b_buf_0.b.a[i]>>1;
p = &(ug->u.a[uId]);
@@ -18553,7 +18702,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
t_ch->b_buf_1.b.n = 0;
get_unitig(ug->g, NULL, bch[1], &convex[1], &bch_occ[1], &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, &(t_ch->b_buf_1));
for (i = 0; i < t_ch->b_buf_1.b.n; ++i)
for (i = 0; i < t_ch->b_buf_1.b.n; ++i)///retrive another side of the bubble
{
uId = t_ch->b_buf_1.b.a[i]>>1;
p = &(ug->u.a[uId]);
@@ -18564,7 +18713,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
t_ch->ir_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
}
}
///generate read-to-read overlaps
chain_origin_trans_uid_s_bubble(&(t_ch->b_buf_0), &(t_ch->b_buf_1),
v0, b->S.a[0]^1, ug, cov);
return;
@@ -31963,6 +32112,10 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
// flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL;
flat_soma_v(sg, sources, ruIndex);
**/
if(!(ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt))) {
// output_unitig_graph(sg, coverage_cut, "pre_clean", sources, ruIndex, max_hang_length, mini_overlap_length);
hic_clean_adv(sg, &uopt);
}
output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources,
asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length);
gfa_ut.cpp
+16 -5
View File
@@ -13127,15 +13127,17 @@ asg64_v *b0, asg64_v *b1, double cutoff)
for (z = 0; z < zn && b0->a[z] == b1->a[z]; z++); ///z: first raw unitig that is different between two paths
assert(z > 0); w0 = w1 = (uint64_t)-1;
if(z < b0->n) get_integer_seq_ovlps(uidx, b0->a, b0->n, z - 1, 0, NULL, &w0);
else return 0;///b0 is contained
if(z < b1->n) get_integer_seq_ovlps(uidx, b1->a, b1->n, z - 1, 0, NULL, &w1);
// if(is_debug) {
else continue;///b1 is contained
// if(((v>>1) == 8722) || ((v>>1) == 56768)) {
// prt_sub_integer_path(b0, it, w0, w1, zn, z);
// prt_sub_integer_path(b1, it, w0, w1, zn, z);
// }
if(w0 == (uint64_t)-1) w0 = 0;
if(w1 == (uint64_t)-1) w1 = 0;
if(b0->n == zn) return 0;///b0 is shorter
if(b1->n == zn) continue;///b1 is shorter
// if(b0->n == zn) return 0;///b0 is shorter
// if(b1->n == zn) continue;///b1 is shorter
if((min_w0 == (uint64_t)-1) || (z == zn) || (min_w0 > w0) || (min_w0 == w0 && min_w1 < w1)) {
min_w0 = w0; min_w1 = w1;
}
@@ -13163,9 +13165,17 @@ uint64_t *ridx, asg64_v *res)
v = int_a[k];
if((!f[v])&&(!f[v^1])) continue;
if((pi != (uint64_t)-1) && (f[v^1])) {
// fprintf(stderr, "+[M::%s::] utg%.6dl(%c), f[v^1]::%u ,v^1::%lu\n",
// __func__, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v^1], v^1);
// if(((v>>1) == 8722) || ((v>>1) == 56768)) {
// fprintf(stderr, "+[M::%s::ii[%lu, %lu)] utg%.6dl(%c), f[v^1]::%u, v^1::%lu, putg%.6dl(%c), f[pv]::%u, pv::%lu\n",
// __func__, s, e, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v^1], v^1,
// (int32_t)(int_a[pi]>>1)+1, "+-"[int_a[pi]&1], f[int_a[pi]], int_a[pi]);
// }
if(is_best_path(uidx, ng, int_idx, int_a, s, e, k, v^1, ridx_a, ridx, b0, b1, 0.51)) {
// if(((v>>1) == 8722) || ((v>>1) == 56768)) {
// fprintf(stderr, "-[M::%s::ii[%lu, %lu)] utg%.6dl(%c), f[v^1]::%u, v^1::%lu, putg%.6dl(%c), f[pv]::%u, pv::%lu\n",
// __func__, s, e, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v^1], v^1,
// (int32_t)(int_a[pi]>>1)+1, "+-"[int_a[pi]&1], f[int_a[pi]], int_a[pi]);
// }
// fprintf(stderr, "[M::%s::] utg%.6dl(%c)->utg%.6dl(%c)\n", __func__,
// (int32_t)(int_a[pi]>>1)+1, "+-"[int_a[pi]&1], (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1]);
pz = (res->n > res_n)? &(res->a[res->n-1]):(NULL);
@@ -16361,4 +16371,5 @@ ul_renew_t *ropt)
ma_hit_contained_advance((*(ropt->src)), (*(ropt->n_read)), (*(ropt->cov)), ropt->ruIndex, ropt->max_hang, ropt->mini_ovlp);
post_rescue(uopt, (*(ropt->sg)), (*(ropt->src)), (*(ropt->r_src)), ropt->ruIndex, ropt->b_mask_t, 0);
// print_raw_uls_aln(uidx, asm_opt.output_file_name);
// exit(0);
}
gfa_ut.h
+3
View File
@@ -1,6 +1,7 @@
#ifndef __GFA_UT__
#define __GFA_UT__
#include "Overlaps.h"
#include "hic.h"
typedef struct {
asg_t *g;
@@ -31,5 +32,7 @@ void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd);
asg_t *gen_ng(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, uint64_t scaffold_len);
void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover);
// 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);
#endif
hic.cpp
+189 -12
View File
@@ -15,6 +15,7 @@
#include "kseq.h" // FASTA/Q parser
#include "kdq.h"
#include "horder.h"
#include "gfa_ut.h"
KSEQ_INIT(gzFile, gzread)
KDQ_INIT(uint64_t)
@@ -16481,15 +16482,14 @@ void optimize_u_trans(kv_u_trans_t *ovlp, kvec_pe_hit* hits, ha_ug_index* idx)
x = &(ovlp->a[i]);
if(x->qn > x->tn) continue;
if(x->f != RC_2 || x->del) continue;
occ = get_oe_occ(x->qn, x->tn, hits, idx) + get_oe_occ(x->tn, x->qn, hits, idx);
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])))));
if(p->nw < 0) fprintf(stderr, "ERROR-nw\n");
if(p->nw == 0) p->nw = x->nw*0.005;
if(p->nw == 0) {
k_trans.n--;
}
else {
} else {
kv_pushp(u_trans_t, k_trans, &p);
(*p) = k_trans.a[k_trans.n-2];
p->qn = k_trans.a[k_trans.n-2].tn; p->qs = k_trans.a[k_trans.n-2].ts; p->qe = k_trans.a[k_trans.n-2].te;
@@ -16540,6 +16540,179 @@ ha_ug_index* idx, uint64_t step, uint64_t total)
}
void prt_hits_noid(ha_ug_index* idx, ma_ug_t* ug, kvec_pe_hit* hits, FILE *fn)
{
uint64_t k, shif = 64 - idx->uID_bits;
char dir[2] = {'+', '-'};
for (k = 0; k < hits->a.n; ++k) {
fprintf(fn, "r-%lu-th\t%c\trs-utg%.6d%c\t%lu\t%c\tre-utg%.6d%c\t%lu\n",
hits->a.a[k].id,
dir[hits->a.a[k].s>>63], (int)((hits->a.a[k].s<<1)>>shif)+1,
"lc"[ug->u.a[((hits->a.a[k].s<<1)>>shif)].circ], hits->a.a[k].s&idx->pos_mode,
dir[hits->a.a[k].e>>63], (int)((hits->a.a[k].e<<1)>>shif)+1,
"lc"[ug->u.a[((hits->a.a[k].e<<1)>>shif)].circ], hits->a.a[k].e&idx->pos_mode);
}
}
void prt_utg_trans(kv_u_trans_t *ta, ma_ug_t* ug, FILE *fn)
{
uint32_t i;
u_trans_t *p = NULL;
for (i = 0; i < ta->n; i++) {
p = &(ta->a[i]);
fprintf(fn, "utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\tw(%f)\tf(%u)\n",
p->qn+1, "lc"[ug->u.a[p->qn].circ], ug->u.a[p->qn].len, p->qs, p->qe, "+-"[p->rev],
p->tn+1, "lc"[ug->u.a[p->tn].circ], ug->u.a[p->tn].len, p->ts, p->te, p->nw, p->f);
}
}
void prt_kv_u_trans(kv_u_trans_t *ta, hc_links* lk, int8_t *s, FILE *fn)
{
uint32_t i;
u_trans_t *p = NULL;
hc_edge *e = NULL;
for (i = 0; i < ta->n; i++) {
p = &(ta->a[i]);
e = get_hc_edge(lk, p->qn, p->tn, 0);
fprintf(fn, "s-utg%.6ul\tS(%d)\td-utg%.6ul\tS(%d)\trev(%u)\td(%lld)\ttw(%f)\n",
p->qn+1, s[p->qn], p->tn+1, s[p->tn], p->rev,
(e == NULL || e->dis == (uint64_t)-1)? -1 : (long long)(e->dis>>3), p->nw);
}
}
void prt_bubble_gfa_adv(FILE *fp, bubble_type *bub, const char* utg_pre, const char* bub_pre, const char* chain_pre)
{
uint32_t i, k, m, *a, n, beg, sink, x; ma_utg_t *p; uint64_t occ;
ma_ug_t *b_ug = bub->b_ug; char name[32], bname[32]; uint8_t *f; CALLOC(f, bub->ug->u.n);
for (i = 0; i < b_ug->u.n; i++) {
p = &b_ug->u.a[i];
if(p->n == 0) continue;
for (k = occ = 0; k < p->n; k++){
x = p->a[k]>>33;
get_bubbles(bub, x, &beg, &sink, &a, &n, NULL);
for (m = 0; m < n; m++) {
occ += bub->ug->u.a[a[m]>>1].n; f[a[m]>>1] = 1;
}
if(beg != (uint32_t)-1 && f[beg>>1] == 0) {
occ += bub->ug->u.a[beg>>1].n; f[beg>>1] = 1;
}
if(sink != (uint32_t)-1 && f[sink>>1] == 0) {
occ += bub->ug->u.a[sink>>1].n; f[sink>>1] = 1;
}
}
sprintf(name, "%s%.6d%c", chain_pre, i + 1, "lc"[p->circ]);
fprintf(fp, "S\t%s\t*\tLN:i:%lu\n", name, occ);
for (k = 0; k < p->n; k++) {
x = p->a[k]>>33;
sprintf(bname, "%s%.6d", bub_pre, x + 1);
fprintf(fp, "B\t%s\t%c\tcid:i:%s\tsm:%c\n", bname, "+-"[(p->a[k]>>32)&1], name, "01"[x<bub->f_bub]);
get_bubbles(bub, x, &beg, &sink, &a, &n, NULL);
if(beg != (uint32_t)-1) {
fprintf(fp, "U\t%s%.6d%c\t%c\tcid:i:%s\tbid:b:%s\thom:%c\n",
utg_pre, (beg>>1)+1, "lc"[bub->ug->u.a[(beg>>1)].circ], "+-"[beg&1], name, bname, "10"[IF_HOM((beg>>1), *bub)]);
}
if(sink != (uint32_t)-1) {
fprintf(fp, "U\t%s%.6d%c\t%c\tcid:i:%s\tbid:s:%s\thom:%c\n",
utg_pre, (sink>>1)+1, "lc"[bub->ug->u.a[(sink>>1)].circ], "+-"[sink&1], name, bname, "10"[IF_HOM((sink>>1), *bub)]);
}
for (m = 0; m < n; m++) {
occ += bub->ug->u.a[a[m]>>1].n;
fprintf(fp, "U\t%s%.6d%c\t%c\tcid:i:%s\tbid:c:%s\thom:%c\n",
utg_pre, (a[m]>>1)+1, "lc"[bub->ug->u.a[(a[m]>>1)].circ], "+-"[a[m]&1], name, bname, "10"[IF_HOM((a[m]>>1), *bub)]);
}
}
}
asg_arc_t* au = NULL;
uint32_t nu, u, v, j;
for (i = 0; i < b_ug->u.n; ++i) {
if(b_ug->u.a[i].m == 0) continue;
if(b_ug->u.a[i].circ)
{
fprintf(fp, "L\t%s%.6dc\t+\t%s%.6dc\t+\t%dM\tL1:i:%d\n",
chain_pre, i+1, chain_pre, i+1, 0, 0);
fprintf(fp, "L\t%s%.6dc\t-\t%s%.6dc\t-\t%dM\tL1:i:%d\n",
chain_pre, i+1, chain_pre, i+1, 0, 0);
}
u = i<<1;
au = asg_arc_a(b_ug->g, u);
nu = asg_arc_n(b_ug->g, u);
for (j = 0; j < nu; j++)
{
if(au[j].del) continue;
v = au[j].v;
fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n",
chain_pre, (u>>1)+1, "lc"[b_ug->u.a[u>>1].circ], "+-"[u&1],
chain_pre, (v>>1)+1, "lc"[b_ug->u.a[v>>1].circ], "+-"[v&1], 0, 0);
}
u = (i<<1) + 1;
au = asg_arc_a(b_ug->g, u);
nu = asg_arc_n(b_ug->g, u);
for (j = 0; j < nu; j++)
{
if(au[j].del) continue;
v = au[j].v;
fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n",
chain_pre, (u>>1)+1, "lc"[b_ug->u.a[u>>1].circ], "+-"[u&1],
chain_pre, (v>>1)+1, "lc"[b_ug->u.a[v>>1].circ], "+-"[v&1], 0, 0);
}
}
for (i = 0; i < bub->ug->u.n; i++) {
if(f[i]) continue;
fprintf(fp, "U\t%s%.6d%c\t+\tcid:i:*\tbid:c:*\thom:%c\n",
utg_pre, i+1, "lc"[bub->ug->u.a[i].circ], "10"[IF_HOM(i, *bub)]);
}
free(f);
}
void prt_debug_hic(const char* o_n, ma_ug_t* ug, ha_ug_index* idx, ug_opt_t *opt, kvec_pe_hit* hits,
kv_u_trans_t *utg_trans, kv_u_trans_t *p_arcs, hc_links* lk, int8_t *s, bubble_type *bub)
{
char* gfa_name = (char*)malloc(strlen(o_n)+100); FILE *fn = NULL;
sprintf(gfa_name, "%s.hic.dbg", o_n);
print_debug_gfa(idx->read_g, idx->ug, opt->coverage_cut, gfa_name, opt->sources, opt->ruIndex,
opt->max_hang, opt->min_ovlp, 0, 0, 0);
if(hits) {
sprintf(gfa_name, "%s.hic.hits.log", o_n); fn = fopen(gfa_name, "w");
prt_hits_noid(idx, ug, hits, fn);
fclose(fn);
}
if(utg_trans) {
sprintf(gfa_name, "%s.hic.utg.trans.log", o_n); fn = fopen(gfa_name, "w");
prt_utg_trans(utg_trans, ug, fn);
fclose(fn);
}
if(p_arcs) {
sprintf(gfa_name, "%s.hic.parcs.log", o_n); fn = fopen(gfa_name, "w");
prt_kv_u_trans(p_arcs, lk, s, fn);
fclose(fn);
}
if(bub) {
sprintf(gfa_name, "%s.bub.noseq.gfa", o_n); fn = fopen(gfa_name, "w");
prt_bubble_gfa_adv(fn, bub, "utg", "btg", "ctg");
fclose(fn);
}
free(gfa_name);
fprintf(stderr, "[M::%s::] done\n", __func__);
}
int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_opt_t *opt, kvec_pe_hit **rhits)
{
double index_time = yak_realtime();
@@ -16608,11 +16781,15 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
// if(bub.round_id == 0) init_phase(idx, &k_trans, &bub, s);
// update_trans_g(idx, &k_trans, &bub);
/*******************************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, (asm_opt.ar)?(&bub):(NULL), &(idx->t_ch->k_trans), 0,
// (((bub.round_id+1) == bub.n_round)?1: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);
(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)**/);
// 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));
// }
/*******************************for debug************************************/
label_unitigs_sm(s->s, NULL, idx->ug);
@@ -16642,15 +16819,15 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
// horder_t *ho = init_horder_t(&sl.hits, idx->uID_bits, idx->pos_mode, idx->read_g, idx->ug, &bub, &(idx->t_ch->k_trans), opt, 3);
///print_hc_links(&link, 0, &hap);
// print_hc_links(&link, 0, &hap);
// print_hits_simp(idx, &sl.hits);
// print_kv_u_trans_t(&(idx->t_ch->k_trans));
// print_kv_u_trans(&k_trans, &link, s->s);
// print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, NULL/**idx->link**/, idx);
// print_hits(idx, &sl.hits, fn1, fn2);
///print_debug_bubble_graph(&bub, idx->ug, asm_opt.output_file_name);
// print_debug_bubble_graph(&bub, idx->ug, asm_opt.output_file_name);
// print_bubble_chain(&bub);
// destory_contig_partition(&hap);
// destory_horder_t(&ho);
kv_destroy(sl.hits.a); kv_destroy(sl.hits.idx); kv_destroy(sl.hits.occ);
rcut.cpp
+598 -2
View File
@@ -1,6 +1,8 @@
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include "assert.h"
#include "rcut.h"
#include "Purge_Dups.h"
#include "Correct.h"
@@ -94,6 +96,28 @@ typedef struct{
bits_p *vis;
}mc_bp_t;
typedef struct{
kvec_t(uint32_t) nn;
kvec_t(uint64_t) ng;
} nn_clus_t;
typedef struct{
bits_p vis;
t_w_t w;
uint32_t off, occ;
}clus_flip_aux;
typedef struct{
nn_clus_t cc;
bubble_type* bub;
const mc_opt_t *opt;
mc_g_t *mg;
uint8_t *lock, lock_max, dbg;
uint32_t n, n_thread;
clus_flip_aux *aux;
mc_svaux_t *baux;
} mc_clus_t;
typedef struct {
uint64_t x; // RNG
uint32_t cc_off, cc_size;
@@ -1685,6 +1709,7 @@ static t_w_t mc_optimize_local(const mc_opt_t *opt, const mc_match_t *ma, mc_sva
if (b->z[k].z[0] == b->z[k].z[1]) continue;
s = b->z[k].z[0] > b->z[k].z[1]? -1 : 1;
if (b->s[k] != s) {
// fprintf(stderr, "utg%.6dl, s[k]::%d, s::%d\n", (int32_t)(k)+1, b->s[k], s);
mc_set_spin(ma, b, k, s);///no need to change the score of k itself
++n_flip;
}
@@ -2209,6 +2234,517 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3
return n_iter;
}
#define clus_a(b, id) (((b)).cc.nn.a+(((b)).cc.ng.a[(id)]>>32))
#define clus_n(b, id) (((uint32_t)((((b)).cc.ng.a[(id)]))))
void reorder_bub(const mc_match_t *ma, uint32_t *a, uint32_t a_n, uint8_t *ff, uint8_t lf, uint8_t rf, uint8_t cf,
uint8_t cuf, double *sc_l, double *sc_r, double *sc_m, uint32_t *res)
{
uint32_t k, o, j, n, t, z; mc_edge_t *e;
double mm_l, mm_r, m_inner, mm; int64_t mmlk, mmrk, rev, m_inner_k, mmk;
for (k = 0; k < a_n; k++) ff[a[k]] = cf;
m_inner = 0; m_inner_k = -1; mm_l = mm_r = 0; mmlk = mmrk = -1;
for (k = 0, rev = -1; k < a_n; k++) {
sc_l[k] = sc_r[k] = sc_m[k] = 0;
o = ma->idx.a[a[k]] >> 32; n = (uint32_t)ma->idx.a[a[k]];
for (j = 0; j < n; ++j) {
e = &ma->ma.a[o + j];
t = ma_y(*e);
if(ff[t] == lf) sc_l[k] += fabs(e->w);
if(ff[t] == rf) sc_r[k] += fabs(e->w);
if(ff[t] == cf) sc_m[k] += fabs(e->w);
}
if(sc_l[k] > 0) {
if(sc_l[k] > mm_l) {
mm_l = sc_l[k]; mmlk = k;
}
} else if(sc_r[k] > 0) {
if(sc_r[k] > mm_r) {
mm_r = sc_r[k]; mmrk = k;
}
} else if(sc_m[k] > 0 && sc_m[k] > m_inner) {
m_inner = sc_m[k]; m_inner_k = k;
}
}
if(mmlk == -1 && mmrk == -1 && m_inner_k == -1) return;
if(mmlk != -1) {
mm = mm_l; mmk = mmlk; rev = 0;
} else if(mmrk != -1) {
mm = mm_r; mmk = mmrk; rev = 1;
} else {
mm = m_inner; mmk = m_inner_k; rev = 0;
}
double *sc = (!rev)?sc_l:sc_r;
// for (k = 0; k < a_n; k++) {
// // sc_m[k] = 0;
// fprintf(stderr, "+[M::%s] a_n::%u, a[%u]::%u\n", __func__, a_n, k, a[k]);
// }
z = 0; res[z++] = a[mmk]; ff[a[mmk]] = cuf;
for (; z < a_n; ) {
for (k = 0, mm = 0, mmk = -1; k < a_n; k++) {
if(ff[a[k]] == cuf) continue;
o = ma->idx.a[a[k]] >> 32; n = (uint32_t)ma->idx.a[a[k]];
for (j = 0; j < n; ++j) {
e = &ma->ma.a[o + j];
t = ma_y(*e);
if(ff[t] == cuf) sc[k] += fabs(e->w);
// if(ff[t] == cuf) sc_m[k] += fabs(e->w);
}
if(sc[k] > 0) {
if(sc[k] > mm) {
mm = sc[k]; mmk = k;
}
}
}
if(mmk == -1) break;
// fprintf(stderr, "-[M::%s] mmk::%ld, a[%ld]::%u, z::%u\n", __func__, mmk, mmk, a[mmk], z);
res[z++] = a[mmk]; ff[a[mmk]] = cuf;
}
// fprintf(stderr, "[M::%s] a_n::%u, z::%u\n", __func__, a_n, k, z);
if(z < a_n) {
for (k = 0; k < a_n; k++) {
if(ff[a[k]] == cuf) continue;
res[z++] = a[k]; ff[a[k]] = cuf;
}
}
assert(z == a_n);
if(!rev) {
for (k = 0; k < a_n; k++) {
a[k] = res[k];
// fprintf(stderr, "[M::%s] a_n::%u, a[k]::%u, res[k]::%u\n", __func__, a_n, a[k], res[k]);
ff[a[k]] = lf;
}
} else {
for (k = 0; k < a_n; k++) {
a[k] = res[a_n-k-1]; ff[a[k]] = lf;
}
}
}
void prt_bub(uint32_t *a, uint32_t a_n, const char *cmd)
{
uint32_t k;
fprintf(stderr, "%s\n", cmd);
for (k = 0; k < a_n; k++) {
fprintf(stderr, "utg%.6dl\t", (int32_t)a[k]+1);
}
fprintf(stderr, "\n");
}
void renew_mc_clus_t(mc_clus_t *bc, uint32_t *a, uint32_t a_n)
{
if(!bc) return;
// fprintf(stderr, "[M::%s] a_n::%u\n", __func__, a_n);
uint32_t k, i, *ba, bn, m, cocc, iin, bub_occ = 0, bbn = 0; uint64_t *p; ma_utg_t *u = NULL;
kvec_t(double) sc_l; kvec_t(double) sc_r; kvec_t(double) sc_m; kvec_t(uint32_t) tmp;
kv_init(sc_l); kv_init(sc_r); kv_init(sc_m); kv_init(tmp);
bc->cc.ng.n = bc->cc.nn.n = 0;
memset(bc->lock, 0, sizeof((*(bc->lock)))*bc->n);
for (k = 0; k < a_n; k++) bc->lock[a[k]] = 1;
kv_resize(uint32_t, bc->cc.nn, a_n);
// for (i = 0; i < bc->bub->chain_weight.n; i++) {
// if(bc->bub->chain_weight.a[i].del) continue;
// u = &(bc->bub->b_ug->u.a[bc->bub->chain_weight.a[i].id]);///list of bubbles
for (i = 0; i < bc->bub->b_ug->u.n; i++) {
u = &(bc->bub->b_ug->u.a[i]);
if(u->n == 0) continue;
bub_occ += u->n;
// fprintf(stderr, "[M::%s] i::%u, u->n::%u\n", __func__, i, (uint32_t)u->n);
for (k = cocc = 0, iin = bc->cc.nn.n; k < u->n; k++) {
get_bubbles(bc->bub, u->a[k]>>33, NULL, NULL, &ba, &bn, NULL);
kv_pushp(uint64_t, bc->cc.ng, &p); bbn += bn;
*p = bc->cc.nn.n;///a bubble
for (m = 0; m < bn; m++) {
if(!(bc->lock[ba[m]>>1])) continue;
kv_push(uint32_t, bc->cc.nn, (ba[m]>>1));
bc->lock[ba[m]>>1] = 2;
}
if(bc->cc.nn.n <= (*p)) {///no node in this bubble
bc->cc.ng.n--;
continue;
}
*p <<= 32; *p |= (bc->cc.nn.n-((*p)>>32)); cocc++;
}
//split chains
if(cocc > 0) {//cocc: # of bubbles in this chain
for (k = bc->cc.ng.n - cocc; k < bc->cc.ng.n; k++) {
kv_resize(double, sc_l, clus_n((*bc), k));
kv_resize(double, sc_r, clus_n((*bc), k));
kv_resize(double, sc_m, clus_n((*bc), k));
kv_resize(uint32_t, tmp, clus_n((*bc), k));
// prt_bub(clus_a((*bc), k), clus_n((*bc), k), "-0-");
reorder_bub(bc->mg->e, clus_a((*bc), k), clus_n((*bc), k), bc->lock, 3, 2, 4, 5,
sc_l.a, sc_r.a, sc_m.a, tmp.a);
// prt_bub(clus_a((*bc), k), clus_n((*bc), k), "-1-");
}
for (k = iin; k < bc->cc.nn.n; k++) bc->lock[bc->cc.nn.a[k]] = 1;//reset
kv_pushp(uint64_t, bc->cc.ng, &p); *p = (uint64_t)-1;
kv_push(uint32_t, bc->cc.nn, ((uint32_t)-1)); ///split
}
}
for (k = 0; k < a_n; k++) bc->lock[a[k]] = 0;
kv_destroy(sc_l); kv_destroy(sc_r); kv_destroy(sc_m); kv_destroy(tmp);
// fprintf(stderr, "[M::%s] a_n::%u, bc->cc.nn.n::%u, bc->cc.ng.n::%u, bbn::%u, bub_occ::%u, bc->bub->b_ug->u.n::%u\n", __func__,
// a_n, (uint32_t)bc->cc.nn.n, (uint32_t)bc->cc.ng.n, bbn, bub_occ, (uint32_t)bc->bub->b_ug->u.n);
}
void clean_clus_flip_aux(clus_flip_aux *z)
{
z->w = -1; z->occ = z->off = (uint32_t)-1;
memset(z->vis.a, 0, sizeof(*(z->vis.a))*z->vis.n);
}
#define is_set_bits_p(v, i) (((v).a[((i)>>3)]>>(i&7))&1)
#define set_bits_p(v, i) (((v).a[((i)>>3)])|=(((uint8_t)1)<<(i&7)));
t_w_t clus_weight(mc_svaux_t *b_aux, mc_match_t *ma, bits_p *vis, uint32_t uid)
{
mc_edge_t *o = NULL;
uint32_t n, i, t;
t_w_t w = ((t_w_t)(b_aux->s[uid])) * (b_aux->z[uid].z[0] - b_aux->z[uid].z[1]) * 2;
t_w_t w_off = 0;
o = pt_a(*ma, uid);
n = pt_n(*ma, uid);
for (i = 0; i < n; ++i) {
t = ma_y(o[i]);
if(!(is_set_bits_p((*vis), t))) continue;
// if(vis[t] == 0) continue;
if(t == uid) continue;
w_off += (b_aux->s[uid]*b_aux->s[t]*o[i].w);
}
return w - (w_off*4);//2 for self; 4 for both directions
}
void cal_clus_sc0(mc_match_t *ma, mc_svaux_t *b_aux, mc_clus_t* bc, uint8_t *lock, uint8_t lock_max,
uint32_t *a, uint32_t a_n, uint32_t id, clus_flip_aux *r, uint32_t tid)
{
uint32_t i, max_occ = (uint32_t)-1;//, len, mm0 = (uint32_t)-1, mm1 = 0;
bits_p *vis = &(r->vis); t_w_t w = 0, max_w = -1;
memset(vis->a, 0, sizeof(*(vis->a))*vis->n);
// if(bc->dbg) {
// if(a[id] == 487 || a[id] == 47) {
// fprintf(stderr, "[M::%s::] id::%u, a[id]::%u, s::%d\n", __func__, id, a[id], b_aux->s[a[id]]);
// }
// }
for (i = id; i < a_n && a[i] != (uint32_t)-1; i++) {
if(lock[a[i]] >= lock_max) continue;
if(is_set_bits_p((*vis), a[i])) continue;
w += clus_weight(b_aux, ma, vis, a[i]);
// if(bc->dbg) {
// if(a[id] == 487 || a[id] == 47) {
// fprintf(stderr, "[M::%s::id->%u] a[%u]::%u, s::%d, lock::%u, lock_max::%u, is_set::%u, w::%f\n",
// __func__, id, i, a[i], b_aux->s[a[i]], lock[a[i]], lock_max, is_set_bits_p((*vis), a[i]), w);
// }
// }
set_bits_p((*vis), a[i]);
// mm1 = a[i]; if(a[i] < mm0) mm0 = a[i];
///update max_w
if(max_w < w) {
max_w = w; max_occ = i + 1 - id;
}
}
// if(mm0 != (uint32_t)-1 && mm1 != (uint32_t)-1) {
// len = MIN(((mm1>>3)+1), vis->n) - (mm0>>3);
// memset(vis->a+(mm0>>3), 0, sizeof(*(vis->a))*len);
// }
for (i = id; i < a_n && a[i] != (uint32_t)-1; i++) vis->a[a[i]>>3] = 0;///reset
if(max_w <= 0.000001 || max_occ == (uint32_t)-1) return;
if((max_w > r->w) || (max_w == r->w && id < r->off)) {
r->w = max_w; r->off = id; r->occ = max_occ;
}
// if(bc->dbg) {
// if(a[id] == 487 || a[id] == 47) {
// fprintf(stderr, "[M::%s::] id::%u, a[id]::%u, w::%f, off::%u, occ::%u\n", __func__, id, a[id], r->w, r->off, r->occ);
// }
// }
}
static void worker_cal_clus_sc(void *data, long i, int tid) // callback for kt_for()
{
mc_clus_t *bc = (mc_clus_t *)data;
uint32_t *a = bc->cc.nn.a, a_n = bc->cc.nn.n;
if(a[i] == (uint32_t)-1) return;
cal_clus_sc0(bc->mg->e, bc->baux, bc, bc->lock, bc->lock_max, a, a_n, i, &(bc->aux[tid]), tid);
}
uint32_t gen_best_clus(mc_clus_t *bc, uint32_t *off, uint32_t *occ, double *rw)
{
uint32_t i; (*off) = (*occ) = (uint32_t)-1; (*rw) = -1;
for (i = 0; i < bc->n_thread; i++) clean_clus_flip_aux(&(bc->aux[i]));
kt_for(bc->n_thread, worker_cal_clus_sc, bc, bc->cc.nn.n);
for (i = 0; i < bc->n_thread; i++) {
if(bc->aux[i].off == (uint32_t)-1) continue;
if(bc->aux[i].w < 0) continue;
if((bc->aux[i].w > (*rw)) || (bc->aux[i].w == (*rw) && bc->aux[i].off < (*off))) {
(*off) = bc->aux[i].off; (*occ) = bc->aux[i].occ; (*rw) = bc->aux[i].w;
}
}
if((*off) != (uint32_t)-1) return 1;
return 0;
}
double flip_chain(const mc_match_t *ma, mc_svaux_t *b, uint32_t *a, uint32_t a_n, uint32_t off, uint32_t occ,
bits_p *vis, uint8_t *lock, uint8_t lock_max)
{
uint32_t k, kn = off + occ;
for (k = off; k < kn; k++) {
vis->a[a[k]>>3] = 0;
// if(occ == 4 && a[off] == 11279) dbg = 1;
}
for (k = off; k < kn; k++) {
if(lock[a[k]] >= lock_max) continue;
if(is_set_bits_p((*vis), a[k])) continue;
set_bits_p((*vis), a[k]); lock[a[k]]++;
// if(dbg) fprintf(stderr, "[M::%s::] a[%u]::%u, s::%d\n", __func__, k, a[k], b->s[a[k]]);
mc_set_spin(ma, b, a[k], -b->s[a[k]]);
}
return mc_score(ma, b);
}
void test_flip_sc(const mc_match_t *ma, mc_svaux_t *b, uint32_t *a, uint32_t a_n, bits_p *vis)
{
mc_edge_t *o = NULL; t_w_t w = 0, w_off = 0;
t_w_t sc_new = mc_score(ma, b), sc;
uint32_t n, i, t, k, uid, z;
for (k = 0; k < a_n; k++) mc_set_spin(ma, b, a[k], -b->s[a[k]]);
for (k = 0; k < a_n; k++) {
uid = a[k];
fprintf(stderr, "+[M::%s::] uid::%u, z[0]::%f, z[1]::%f\n", __func__, uid, b->z[uid].z[0], b->z[uid].z[1]);
w += ((t_w_t)(b->s[uid])) * (b->z[uid].z[0] - b->z[uid].z[1]) * 2;
o = pt_a(*ma, uid);
n = pt_n(*ma, uid);
for (i = 0; i < n; ++i) {
t = ma_y(o[i]);
for (z = 0; z < k; z++) {
if(a[z] == t) break;
}
if(z >= k) continue;
// if(!(is_set_bits_p((*vis), t))) continue;
// if(vis[t] == 0) continue;
if(t == uid) continue;
fprintf(stderr, "+[M::%s::] uid::%u, t::%u, w::%f\n", __func__, uid, t, o[i].w);
w_off += (b->s[uid]*b->s[t]*o[i].w);
}
}
w -= (w_off*4);
sc = mc_score(ma, b);
fprintf(stderr, "+[M::%s::] sc::%f, sc_new::%f, w::%f\n", __func__, sc, sc_new, w);
w = w_off = 0; memset(vis->a, 0, sizeof(*(vis->a))*vis->n);
for (k = 0; k < a_n; k++) {
uid = a[k];
fprintf(stderr, "-[M::%s::] uid::%u, z[0]::%f, z[1]::%f\n", __func__, uid, b->z[uid].z[0], b->z[uid].z[1]);
w += ((t_w_t)(b->s[uid])) * (b->z[uid].z[0] - b->z[uid].z[1]) * 2;
o = pt_a(*ma, uid);
n = pt_n(*ma, uid);
for (i = 0; i < n; ++i) {
t = ma_y(o[i]);
// for (z = 0; z < k; z++) {
// if(a[z] == t) break;
// }
// if(z >= k) continue;
if(!(is_set_bits_p((*vis), t))) continue;
// if(vis[t] == 0) continue;
if(t == uid) continue;
fprintf(stderr, "-[M::%s::] uid::%u, t::%u, w::%f\n", __func__, uid, t, o[i].w);
w_off += (b->s[uid]*b->s[t]*o[i].w);
}
set_bits_p((*vis), uid);
}
w -= (w_off*4);
sc = mc_score(ma, b);
fprintf(stderr, "-[M::%s::] sc::%f, sc_new::%f, w::%f\n", __func__, sc, sc_new, w);
}
t_w_t mc_solve_clus(mc_clus_t *bc)
{
uint32_t off, occ/**, r = 0**/; double rw, sc, sc_opt = mc_score(bc->mg->e, bc->baux);
memset(bc->lock, 0, bc->n*sizeof(*(bc->lock)));
while (gen_best_clus(bc, &off, &occ, &rw)) {
sc = flip_chain(bc->mg->e, bc->baux, bc->cc.nn.a, bc->cc.nn.n, off, occ, &(bc->aux[0].vis), bc->lock, bc->lock_max);
// if(r%10000) {
// fprintf(stderr, "[M::%s::] rw::%f, sc_opt::%f, sc::%f, off::%u, occ::%u, r::%u\n",
// __func__, rw, sc_opt, sc, off, occ, r);
// }
if(sc < sc_opt) {
fprintf(stderr, "\nwrong::[M::%s::] rw::%f, sc_opt::%f, sc::%f, off::%u, occ::%u\n",
__func__, rw, sc_opt, sc, off, occ);
// if(occ == 2) {
// uint32_t k;
// for (k = off; k < off + occ; k++) {
// fprintf(stderr, "[M::%s::] a[%u]::%u, s::%d\n", __func__, k, bc->cc.nn.a[k], bc->baux->s[bc->cc.nn.a[k]]);
// }
// test_flip_sc(bc->mg->e, bc->baux, bc->cc.nn.a+off, occ, &(bc->aux[0].vis));
// }
}
sc_opt = sc; //r++;
}
return sc_opt;
}
t_w_t mc_clus_cc(mc_clus_t *bc)
{
// fprintf(stderr, "+[M::%s::]\n", __func__);
// double index_time = yak_realtime();
// uint32_t r = 1;
t_w_t sc_opt, sc;
// fprintf(stderr, "-[M::%s::]\n", __func__);
mc_reset_z(bc->mg->e, bc->baux);
// fprintf(stderr, "*[M::%s::]\n", __func__);
sc_opt = mc_score(bc->mg->e, bc->baux);
// fprintf(stderr, "[M::%s::] sc_opt: %f\n", __func__, sc_opt);
while (1) {
sc = mc_solve_clus(bc);
// fprintf(stderr, "[M::%s::# round: %u] sc_opt: %f, sc: %f\n", __func__, r, sc_opt, sc);
if(sc <= (sc_opt+0.0000001)) break;
sc_opt = sc; //r++;
}
// fprintf(stderr, "[M::%s::%.3f] ==> round %u\n", __func__, yak_realtime()-index_time, r);
return sc;
}
uint32_t mc_solve_cc_adv(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint32_t cc_off, uint32_t cc_size, mc_clus_t *bc)
{
uint32_t j, k, n_iter = 0, flush = opt->max_iter * 50, n_skip, n_skip_flush = opt->n_perturb/16;
t_w_t sc_opt = -(1<<30), sc;///problem-w
b->cc_off = cc_off, b->cc_size = cc_size;
if (b->cc_size < 2) return 0;
sc_opt = mc_init_spin(mg->e, b);
// print_sc(opt, mg, b, sc_opt, n_iter);
if (b->cc_size == 2) return 0;
for (j = 0; j < b->cc_size; ++j) {///backup s and z in s_opt and z_opt
b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; ///hap status of each unitig
b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; ///z[0]: positive weight; z[1]: positive weight
}
renew_mc_clus_t(bc, b->cc_node, b->cc_size);
// fprintf(stderr, "\ncc_size: %u, cc_off: %u\n", b->cc_size, b->cc_off);
// print_sc(opt, mg->e, b, sc_opt, n_iter);
sc = mc_optimize_local(opt, mg->e, b, &n_iter);
if (sc > sc_opt) {
for (j = 0; j < b->cc_size; ++j) {
b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]];
b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]];
}
sc_opt = sc;
} else {
for (j = 0; j < b->cc_size; ++j) {
b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]];
b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]];
}
}
// print_mc_node(mg->e, b, 36880);
// print_sc(opt, mg, b, sc_opt, n_iter);
// mc_reset_z_debug(mg->e, b);
// print_sc(opt, mg->e, b, sc_opt, n_iter);
// fprintf(stderr, "\ncc_size: %u, cc_off: %u\n", b->cc_size, b->cc_off);
if(bc) {
sc = mc_clus_cc(bc);
if (sc > sc_opt) {
for (j = 0; j < b->cc_size; ++j) {
b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]];
b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]];
}
sc_opt = sc;
} else {
for (j = 0; j < b->cc_size; ++j) {
b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]];
b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]];
}
}
}
for (k = n_skip = 0; k < (uint32_t)opt->n_perturb; ++k) {
if (k&1) mc_perturb(opt, mg->e, b);
else mc_perturb_node(opt, mg->e, b, 3);
sc = mc_optimize_local(opt, mg->e, b, &n_iter);
// if((k%256) == 0) fprintf(stderr, "(%u) sc_opt::%f, sc::%f, flush::%u, n_iter::%u\n", k, sc_opt, sc, flush, n_iter);
if (sc > sc_opt) {
for (j = 0; j < b->cc_size; ++j) {
b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]];
b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]];
}
sc_opt = sc; n_skip = 0;
// print_sc(opt, mg, b, sc_opt, n_iter);
} else {
for (j = 0; j < b->cc_size; ++j) {
b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]];
b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]];
}
n_skip++;
}
if(n_skip >= n_skip_flush && bc) {
sc = mc_clus_cc(bc);
if (sc > sc_opt) {
for (j = 0; j < b->cc_size; ++j) {
b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]];
b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]];
}
sc_opt = sc;
} else {
for (j = 0; j < b->cc_size; ++j) {
b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]];
b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]];
}
}
n_skip = 0;
}
if((n_iter%flush) == 0) {
mc_reset_z(mg->e, b);
sc = mc_score(mg->e, b);
for (j = 0; j < b->cc_size; ++j) {
b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]];
}
sc_opt = sc;
}
}
if(bc) {
sc = mc_clus_cc(bc);
if (sc > sc_opt) {
for (j = 0; j < b->cc_size; ++j) {
b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]];
b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]];
}
sc_opt = sc;
} else {
for (j = 0; j < b->cc_size; ++j) {
b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]];
b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]];
}
}
}
// if(bc) {
// bc->dbg = 1;
// mc_clus_cc(bc);
// bc->dbg = 0;
// }
for (j = 0; j < b->cc_size; ++j)
{
b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]];
b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]];
}
// exit(1);
return n_iter;
}
void reset_mb_g_t_z(mb_g_t *mbg);
uint32_t mb_solve_cc(const mc_opt_t *opt, mb_g_t *mbg, mb_svaux_t *b, uint32_t cc_off, uint32_t cc_size)
@@ -2689,6 +3225,63 @@ void mc_solve_core(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub)
fprintf(stderr, "[M::%s::%.3f] ==> Partition\n", __func__, yak_realtime()-index_time);
}
mc_clus_t *init_mc_clus_t(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub, uint32_t n_thread, mc_svaux_t *b, uint32_t flp_max)
{
if((!bub)) return NULL;
mc_clus_t *p; CALLOC(p, 1);
p->bub = bub; p->opt = opt; p->mg = mg; p->n = bub->ug->g->n_seq;
CALLOC(p->lock, p->n);
if(n_thread > 64) n_thread = 64; p->n_thread = n_thread;
CALLOC(p->aux, p->n_thread);
uint32_t k, ss = (p->n>>3)+(!!(p->n&7));
for (k = 0; k < p->n_thread; k++) {
kv_resize(uint8_t, p->aux[k].vis, ss); p->aux[k].vis.n = ss;
}
p->baux = b; p->lock_max = ((flp_max<=255)?flp_max:255); if(p->lock_max < 1) p->lock_max = 1;
return p;
}
void mc_solve_core_adv(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub)
{
double index_time = yak_realtime();
uint32_t st, i;
mc_svaux_t *b; mc_clus_t *bc;
// mc_bp_t *bp = NULL;
mc_g_cc(mg->e);
b = mc_svaux_init(mg, opt->seed);
bc = init_mc_clus_t(opt, mg, bub, asm_opt.thread_num, b, 16);
// bc = gen_mc_clus_t(mg->e, b, bub, ref, asm_opt.thread_num);
// if(bub) bp = mc_bp_t_init(mg->e, b, bub, asm_opt.thread_num);
/*******************************for debug************************************/
// if(bp) mc_init_spin_all(opt, mg, NULL, b);
// if(bp) mc_solve_bp(bp);
/*******************************for debug************************************/
if(VERBOSE_CUT)
{
fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
}
for (st = 0, i = 1; i <= mg->e->n_seq; ++i) {
if (i == mg->e->n_seq || mg->e->cc[st]>>32 != mg->e->cc[i]>>32) {
mc_solve_cc_adv(opt, mg, b, st, i - st, bc);
st = i;
}
}
if(VERBOSE_CUT)
{
fprintf(stderr, "##############end-[---M::%s::score->%f] ==> Partition\n", __func__, mc_score_all(mg->e, b));
mc_status_all(mg->e, mg->s.a);
}
// if(bp) mc_solve_bp(bp);
///mc_write_info(g, b);
mc_svaux_destroy(b);
// if(bp) destroy_mc_bp_t(&bp);
fprintf(stderr, "[M::%s::%.3f] ==> Partition\n", __func__, yak_realtime()-index_time);
}
void set_p_flag(mc_g_t *mg, uint32_t uID, uint8_t* trio_flag, trans_chain* t_ch, int8_t s)
{
uint32_t i;
@@ -2908,6 +3501,7 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u
{
if(is_dump) {
dump_debug_phasing(MC_NAME, ta, ug, read_g, f_rate, renew_s, s, is_sys, bub, ref);
// bub = NULL;
}
mc_opt_t opt;
@@ -2920,7 +3514,8 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u
mb_solve_core(&opt, mg, ref, is_sys);
///debug_mc_g_t(mg);
// if(renew_s == 0) write_mc_g_t(&opt, mg, MC_NAME);
mc_solve_core(&opt, mg, bub);
// mc_solve_core(&opt, mg, bub);
mc_solve_core_adv(&opt, mg, bub);
if((asm_opt.flag & HA_F_PARTITION) && t_ch)
{
@@ -4005,9 +4600,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)\n", (int32_t)(k)+1, ug->g->seq[k].len);
fprintf(stderr, "utg%.6dl(len::%u), s[k]::%d\n", (int32_t)(k)+1, ug->g->seq[k].len, s[k]);
}