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

update for rescuing missing contigs

Browse Source
This commit is contained in:
chhylp123
2023-04-06 23:14:41 -04:00
parent 503cd0f7cc
commit 5df5892f33
5 changed files with 490 additions and 26 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.h
+1 -1
View File
@@ -5,7 +5,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.19.3-r572"
#define HA_VERSION "0.19.4-r575"
#define VERBOSE 0
Overlaps.cpp
+141 -18
View File
@@ -15208,11 +15208,10 @@ void purge_ovlp_cov_adv(uint32_t id, kv_u_trans_t *ta, asg64_v *b64, ug_rid_cov_
}
}
void trans_sec_cut_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc* src)
void trans_sec_cut_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc* src, ug_rid_cov_t *in)
{
uint64_t k;
asg64_v b64; kv_init(b64);
ug_rid_cov_t *cc = gen_ug_rid_cov_t(ug, sg, src);
uint64_t k; asg64_v b64; kv_init(b64);
ug_rid_cov_t *cc = ((in)?(in):(gen_ug_rid_cov_t(ug, sg, src)));
fprintf(stderr, "+[M::%s]\thom_cov::%lu\thet_cov::%lu\thom_cut::%lu\n", __func__,
cc->hom_cov, cc->het_cov, cc->hom_max);
@@ -15221,7 +15220,7 @@ void trans_sec_cut_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *sg, ma
purge_ovlp_cov_adv(k, ta, &b64, cc, cc->hom_max);
}
destory_ug_rid_cov_t(cc); free(cc); kv_destroy(b64);
if(!in) {destory_ug_rid_cov_t(cc); free(cc);} kv_destroy(b64);
}
static void worker_for_trans_sec_simple_cut(void *data, long i, int tid) // callback for kt_for()
@@ -15251,8 +15250,18 @@ static void worker_for_trans_sec_simple_cut(void *data, long i, int tid) // call
}
}
void gen_ug_rid_cov_t_by_ovlp(kv_u_trans_t *ta, ug_rid_cov_t *cc)
{
uint64_t z, k, id, n; u_trans_t *a;
for (z = 0; z < cc->ug->g->n_seq; z++) {
id = z; a = u_trans_a(*ta, id); n = u_trans_n(*ta, id);
for (k = 0; k < n; k++) {
append_cov_line_ug_rid_cov_t(id, cc->cov.a+cc->idx[id], &(a[k]), cc, ((uint64_t)-1), -1);
}
}
}
void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* src)
void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* src, ug_rid_cov_t *in)
{
u_trans_clean_t sl; uint64_t k, i, l, st, occ; ha_mzl_t *tz;
ha_mzl_v srt_a; kv_u_trans_t *bl; u_trans_t *z;
@@ -15318,7 +15327,7 @@ void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *
free(sl.res); kv_destroy(srt_a);
kt_u_trans_t_idx(ta, ug->g->n_seq);
// dbg_prt_utg_trans(ta, ug, "after");
trans_sec_cut_filter_mmhap_adv(ta, ug, read_g, src);
trans_sec_cut_filter_mmhap_adv(ta, ug, read_g, src, in);
kt_for(sl.n_thread, worker_for_trans_sec_simple_cut, &sl, sl.ta->idx.n);
// CALLOC(sl.srt, sl.n_thread); sl.sec_rate = 0.5;
@@ -15338,6 +15347,7 @@ void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *
st = i;
}
}
// dbg_prt_utg_extra_trans(ta, ug, asm_opt.output_file_name);
}
@@ -15660,9 +15670,9 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, long long gap
// if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(cov->t_ch, output_file_name);
}
dbg_prt_utg_trans(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, "pre");
clean_u_trans_t_idx_filter_mmhap_adv(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, sg, opt->sources);
dbg_prt_utg_trans(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, "after");
// dbg_prt_utg_trans(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, "pre");
clean_u_trans_t_idx_filter_mmhap_adv(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, sg, opt->sources, NULL);
// dbg_prt_utg_trans(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, "after");
// refine_hic_trans_mmhap(opt, &(cov?cov->t_ch->k_trans:t_ch->k_trans), sg, ug);
///for debug
@@ -15723,7 +15733,7 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, long long gap
kv_destroy(d_edges.a);
asg_cleanup(sg);
dbg_prt_trio_mmhap_label(ug, rh, output_file_name);
// dbg_prt_trio_mmhap_label(ug, rh, output_file_name);
output_trio_mmhap(sg, coverage_cut, output_file_name, sources, reverse_sources, tipsLen, tip_drop_ratio,
stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, gap_fuzz, b_mask_t, opt, ug, rh, asm_opt.polyploidy);
@@ -18596,46 +18606,62 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
{
asg_t *g = ug->g;
uint32_t is_first = 1;
// if(trio_flag == MOTHER) {
// print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len, 0, 1, 0);
// exit(1);
// }
redo:
///print_untig((ug), 61955, "i-0:", 0);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-0:", 0);
// fprintf(stderr, "[M::%s] 0\n", __func__);
asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, NULL, 1);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-1:", 0);
// fprintf(stderr, "[M::%s] 1\n", __func__);
magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate);
// fprintf(stderr, "[M::%s] 2\n", __func__);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-2:", 0);
/**********debug**********/
if(just_bubble_pop == 0)
{
cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
}
// fprintf(stderr, "[M::%s] 3\n", __func__);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-3:", 0);
/**********debug**********/
long long pre_cons = get_graph_statistic(g);
long long cur_cons = 0;
while(pre_cons != cur_cons)
{
// fprintf(stderr, "[M::%s] 4\n", __func__);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-4:", 0);
pre_cons = get_graph_statistic(g);
// fprintf(stderr, "[M::%s] 5\n", __func__);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-5:", 0);
///need consider tangles
asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, NULL, 1);
// fprintf(stderr, "[M::%s] 6\n", __func__);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-6:", 0);
/**********debug**********/
if(just_bubble_pop == 0)
{
///need consider tangles
asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, trio_flag, cov, NULL);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-7:", 0);
// fprintf(stderr, "[M::%s] 7\n", __func__);
// if(trio_flag == MOTHER) print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, cov, NULL);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-8:", 0);
// fprintf(stderr, "[M::%s] 8\n", __func__);
asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL, trio_flag);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-9:", 0);
// fprintf(stderr, "[M::%s] 9\n", __func__);
asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL, trio_flag);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-10:", 0);
// fprintf(stderr, "[M::%s] 10\n", __func__);
detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL, cov->is_r_het);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-11:", 0);
// fprintf(stderr, "[M::%s] 11\n", __func__);
///need consider tangles
///note we need both the read graph and the untig graph
@@ -18643,34 +18669,42 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
/**********debug**********/
cur_cons = get_graph_statistic(g);
// fprintf(stderr, "[M::%s] 12\n", __func__);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-12:", 0);
}
if(just_bubble_pop == 0)
{
// fprintf(stderr, "[M::%s] 13\n", __func__);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-13:", 0);
cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-14:", 0);
// fprintf(stderr, "[M::%s] 14\n", __func__);
}
// print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
// fprintf(stderr, "[M::%s] 15\n", __func__);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-15:", 0);
magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate);
// fprintf(stderr, "[M::%s] 16\n", __func__);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-16:", 0);
// print_debug_gfa(read_g, ug, coverage_cut, "resolve_tangles", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len, 0, 0, 0);
// exit(1);
///bug here
resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, trio_flag, drop_ratio);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-17:", 0);
// fprintf(stderr, "[M::%s] 17\n", __func__);
drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-18:", 0);
// fprintf(stderr, "[M::%s] 18\n", __func__);
all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, cov->is_r_het, DOUBLE_CHECK_THRES, asm_opt.trio_flag_occ_thres);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-19:", 0);
// fprintf(stderr, "[M::%s] 19\n", __func__);
// if(trio_flag == MOTHER) print_untig_by_read(ug, "m54329U_190827_173812/30214441/ccs", (uint32_t)-1, NULL, NULL, "bf-16");
if(is_first)
{
is_first = 0;
unitig_arc_del_short_diploid_by_length(ug->g, drop_ratio);
// if(trio_flag == MOTHER) print_untig((ug), 425, "i-20:", 0);
// fprintf(stderr, "[M::%s] 20\n", __func__);
goto redo;
}
@@ -19332,16 +19366,105 @@ void delete_useless_nodes(ma_ug_t **ug)
asg_cleanup(nsg);
}
inline uint32_t is_useful_node(uint32_t flag_occ, uint32_t non_flag_occ, uint32_t drop_occ, uint32_t tot_occ,
float flag_rate, float used_rate, uint32_t min_occ)
{
if((flag_occ > 0) && (flag_occ >= min_occ) && (flag_occ >= ((non_flag_occ+flag_occ+drop_occ)*flag_rate))
&& (drop_occ <= (tot_occ*used_rate))) {
return 1;
}
return 0;
}
void recover_chain_nodes(buf_t *in, ma_ug_t *ug, uint32_t flag, float flag_rate, float used_rate, uint32_t min_occ)
{
ma_utg_t *u = NULL; uint32_t flag_occ, non_flag_occ, drop_occ, rid;
uint32_t tot_flag_occ, tot_non_flag_occ, tot_drop_occ, tot_occ, i, k, z;
tot_flag_occ = tot_non_flag_occ = tot_drop_occ = tot_occ = 0;
for (i = 0; i < in->b.n; i++) {
u = &(ug->u.a[in->b.a[i]>>1]);
flag_occ = non_flag_occ = drop_occ = 0;
for (k = 0; k < u->n; k++) {
rid = u->a[k]>>33;
if(R_INF.trio_flag[rid] == AMBIGU) continue;
else if(R_INF.trio_flag[rid] == DROP) drop_occ++;
else if(R_INF.trio_flag[rid] == flag) flag_occ++;
else if(R_INF.trio_flag[rid] != flag) non_flag_occ++;
}
if(is_useful_node(flag_occ, non_flag_occ, drop_occ, u->n, flag_rate, used_rate, min_occ)) {
ug->g->seq[in->b.a[i]>>1].c = PRIMARY_LABLE;
}
tot_flag_occ += flag_occ;
tot_non_flag_occ += non_flag_occ;
tot_drop_occ += drop_occ;
tot_occ += u->n;
if(is_useful_node(tot_flag_occ, tot_non_flag_occ, tot_drop_occ, tot_occ, flag_rate, used_rate, min_occ)) {
for (z = 0; z <= i; z++) {
if(ug->g->seq[in->b.a[z]>>1].c == ALTER_LABLE) {
u = &(ug->u.a[in->b.a[z]>>1]);
flag_occ = non_flag_occ = drop_occ = 0;
for (k = 0; k < u->n; k++) {
rid = u->a[k]>>33;
if(R_INF.trio_flag[rid] == AMBIGU) continue;
else if(R_INF.trio_flag[rid] == DROP) drop_occ++;
else if(R_INF.trio_flag[rid] == flag) flag_occ++;
else if(R_INF.trio_flag[rid] != flag) non_flag_occ++;
}
if(is_useful_node(flag_occ, non_flag_occ, drop_occ, u->n, flag_rate, used_rate, 0/**min_occ**/)) {
ug->g->seq[in->b.a[z]>>1].c = PRIMARY_LABLE;
}
}
}
}
}
}
void rescue_useless_trio_nodes(ma_ug_t *ug, uint32_t flag, float flag_rate, float used_rate, uint32_t min_occ)
{
asg_t* nsg = ug->g;
uint32_t v, w, n_vtx = nsg->n_seq<<1, i, k, z;
long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen;
buf_t b; memset(&b, 0, sizeof(buf_t));
for (v = 0; v < n_vtx; v++) {
if(nsg->seq[v>>1].del) continue;
if(nsg->seq[v>>1].c != ALTER_LABLE) continue;
///check if beg is the tig end
if(get_real_length(nsg, v^1, NULL) == 1) {
get_real_length(nsg, v^1, &w);
if(get_real_length(nsg, w^1, NULL) == 1) continue;
}
b.b.n = 0;
get_unitig(nsg, ug, v, &w, &nodeLen, &baseLen, &max_stop_nodeLen, &max_stop_baseLen, 1, &b);
recover_chain_nodes(&b, ug, flag, flag_rate, used_rate, min_occ);
if(b.b.n > 1) {
k = b.b.n>>1;
for (i = 0; i < k; i++) {
z = b.b.a[i]; b.b.a[i] = b.b.a[b.b.n-i-1]; b.b.a[b.b.n-i-1] = z;
}
recover_chain_nodes(&b, ug, flag, flag_rate, used_rate, min_occ);
}
}
free(b.b.a);
asg_cleanup(nsg);
}
void delete_useless_trio_nodes(ma_ug_t **ug, asg_t* read_g, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, R_to_U* ruIndex)
ma_hit_t_alloc* sources, R_to_U* ruIndex, uint32_t flag, float flag_rate, float used_rate, uint32_t min_occ)
{
asg_t* nsg = (*ug)->g;
uint32_t v, n_vtx = nsg->n_seq;
uint8_t* primary_flag = (uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t));
if(flag_rate > 0 && used_rate > 0 && min_occ > 0) {
rescue_useless_trio_nodes(*ug, flag, flag_rate, used_rate, min_occ);
}
for (v = 0; v < n_vtx; ++v)
{
if(nsg->seq[v].del) continue;
@@ -19630,7 +19753,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
// fprintf(stderr, "[M::%s] 4\n", __func__);
///delete_useless_nodes(ug);
delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex);
delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex, flag, 0.8, 0.15, 16);
// fprintf(stderr, "[M::%s] 5\n", __func__);
update_hap_label(*ug, read_g);
@@ -19678,7 +19801,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
// fprintf(stderr, "[M::%s] 16\n", __func__);
///delete_useless_nodes(ug);
delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex);
delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex, flag, 0.8, 0.15, 16);
// fprintf(stderr, "[M::%s] 17\n", __func__);
@@ -19688,7 +19811,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 1, 0,
cov, 0, 0);
///delete_useless_nodes(ug);
delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex);
delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex, flag, 0.8, 0.15, 16);
}
// fprintf(stderr, "[M::%s] 18\n", __func__);
Overlaps.h
+2
View File
@@ -1233,5 +1233,7 @@ void destory_ug_rid_cov_t(ug_rid_cov_t *p);
uint32_t append_cov_line_ug_rid_cov_t(uint64_t uid, uint64_t *qcc, u_trans_t *p, ug_rid_cov_t *idx, uint64_t hom_cut, double cut_rate);
uint64_t infer_mmhap_copy(ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc *src, uint8_t *ff, uint64_t uid, uint64_t het_cov, uint64_t n_hap);
uint64_t trans_sec_cut0(kv_u_trans_t *ta, asg64_v *srt, uint32_t id, double sec_rate, uint64_t bd, ma_ug_t *ug);
void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* src, ug_rid_cov_t *in);
void gen_ug_rid_cov_t_by_ovlp(kv_u_trans_t *ta, ug_rid_cov_t *cc);
#endif
README.md
+11 -7
View File
@@ -43,11 +43,8 @@ See [tutorial][tutorial] for more details.
## <a name="intro"></a>Introduction
Hifiasm is a fast haplotype-resolved de novo assembler for PacBio HiFi reads.
It can assemble a human genome in several hours and assemble a ~30Gb California
redwood genome in a few days. Hifiasm emits partially phased assemblies of
quality competitive with the best assemblers. Given parental short reads or
Hi-C data, it produces arguably the best haplotype-resolved assemblies so far.
Hifiasm is a fast haplotype-resolved de novo assembler initially designed for PacBio HiFi reads.
Its latest release could support the telomere-to-telomere assembly by utilizing ultralong Oxford Nanopore reads. Hifiasm produces arguably the best single-sample telomere-to-telomere assemblies combing HiFi, ultralong and Hi-C reads, and it is one of the best haplotype-resolved assemblers for the trio-binning assembly given parental short reads. For a human genome, hifiasm can produce the telomere-to-telomere assembly in one day.
## <a name="why"></a>Why Hifiasm?
@@ -146,11 +143,18 @@ The second command line will run much faster than the first.
### <a name="ul"></a>Ultra-long ONT integration
Hifiasm could integrate ultra-long ONT reads to improve the assembly quality:
Hifiasm could integrate ultra-long ONT reads to produce the telomere-to-telomere assembly:
```sh
hifiasm -o NA12878.asm -t32 --ul ul.fq.gz HiFi-reads.fq.gz
```
Please note that this mode is not stable right now. We have only tested with >=100kb UL reads.
For the single-sample telomere-to-telomere assembly with Hi-C reads:
```sh
hifiasm -o NA12878.asm -t32 --ul ul.fq.gz --h1 read1.fq.gz --h2 read2.fq.gz HiFi-reads.fq.gz
```
For the trio-binning telomere-to-telomere assembly；
```sh
hifiasm -o NA12878.asm -t32 --ul ul.fq.gz -1 pat.yak -2 mat.yak HiFi-reads.fq.gz
```
### <a name="output"></a>Output files
inter.cpp
+335
View File
@@ -10434,6 +10434,202 @@ static void worker_for_trans_ovlp_mmhap_adv(void *data, long i, int tid) // call
s->free_cnt[tid]++;
}
uint32_t tranfor_ovlp(u_trans_t *qovlp, u_trans_t *tovlp, asg_t *g, ul_ov_t *res, uint32_t adjust_rev)
{
int64_t os, oe, s_shift, e_shift, tt, qs, qe, ts, te;
os = MAX(qovlp->ts, tovlp->qs);
oe = MIN(qovlp->te, tovlp->qe);
if(oe <= os) return 0;
///[os, oe) -> qovlp->t*
s_shift = get_offset_adjust(os-qovlp->ts, qovlp->te-qovlp->ts, qovlp->qe-qovlp->qs);
e_shift = get_offset_adjust(qovlp->te-oe, qovlp->te-qovlp->ts, qovlp->qe-qovlp->qs);
if(qovlp->rev) {
tt = s_shift; s_shift = e_shift; e_shift = tt;
}
qs = qovlp->qs+s_shift; qe = ((int64_t)qovlp->qe)-e_shift;
if(qs >= qe) return 0;
///[os, oe) -> tovlp->q*
s_shift = get_offset_adjust(os-tovlp->qs, tovlp->qe-tovlp->qs, tovlp->te-tovlp->ts);
e_shift = get_offset_adjust(tovlp->qe-oe, tovlp->qe-tovlp->qs, tovlp->te-tovlp->ts);
if(tovlp->rev) {
tt = s_shift; s_shift = e_shift; e_shift = tt;
}
ts = tovlp->ts+s_shift; te = ((int64_t)tovlp->te)-e_shift;
if(ts >= te) return 0;
memset(res, 0, sizeof(*res));
res->qn = qovlp->qn; res->qs = qs; res->qe = qe;
res->tn = tovlp->tn; res->ts = ts; res->te = te;
res->rev = ((qovlp->rev == tovlp->rev)?0:1);
if(adjust_rev && res->rev) {///for linear chaining
res->ts = g->seq[res->tn].len - te;
res->te = g->seq[res->tn].len - ts;
}
return 1;
}
uint32_t rescue_adject_ovlp(asg_t *g, uint32_t id, kv_u_trans_t *ta, kv_ul_ov_t *out, st_mt_t *buf)
{
u_trans_t *a, *b; ul_ov_t rr; uint64_t a_n, b_n, k, l, i, z, m;
a = u_trans_a((*ta), id); a_n = u_trans_n((*ta), id);
for (i = out->n = buf->n = 0; i < a_n; i++) {
b = u_trans_a((*ta), a[i].tn); b_n = u_trans_n((*ta), a[i].tn);
z = a[i].tn; z <<= 32; kv_push(uint64_t, *buf, z);
for (k = 0; k < b_n; k++) {
if(b[k].tn == id) continue;
if(!tranfor_ovlp(&(a[i]), &(b[k]), g, &rr, 1)) continue;
z = rr.tn; z <<= 32; z |= out->n; z |= ((uint64_t)0x80000000);
rr.tn <<= 1; rr.tn |= rr.rev; kv_push(ul_ov_t, *out, rr);
}
}
if(out->n == 0) return 1;
radix_sort_gfa64(buf->a, buf->a + buf->n);
for (k = 1, l = m = 0; k <= buf->n; k++) {
if(k == buf->n || (buf->a[l]>>32)!=(buf->a[k]>>32)) {
if((k - l > 1) && (!(buf->a[l]&((uint64_t)0x80000000)))) {///overlap within bck
for (z = l; z < k; z++) {
if(buf->a[z]&((uint64_t)0x80000000)) {
out->a[(uint32_t)(buf->a[z]-((uint64_t)0x80000000))].tn = (uint32_t)-1;
m++;
}
}
}
l = k;
}
}
if(m) {
for (k = m = 0; k < out->n; k++) {
if(out->a[k].tn == (uint32_t)-1) continue;
out->a[m++] = out->a[k];
}
out->n = m;
}
if(out->n == 0) return 1;
radix_sort_ul_ov_srt_tn(out->a, out->a+out->n);
for (k = 0; k < out->n; k++) out->a[k].tn >>= 1;
return 0;
}
/**
uint64_t gen_trans_chain_mmhap(ug_trans_t *s, uint64_t rid, ha_ovec_buf_t *b, kv_ul_ov_t *bl, char *seq, uint64_t len,
double err, double bw)
{
uint64_t cnt = ((s->idx_n.a[rid+1]-s->idx_n.a[rid])), ol_h = 0, pass_aln = 0;
uint32_t high_occ = asm_opt.polyploidy + 1; overlap_region *aux_o = NULL;
///note: high_occ is different
ug_map_lchain(b->abl, rid, seq, len, s->w, s->k, &(s->udb), &b->olist, &b->clist, bw, bw,
s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3,
s->is_HPC, s->idx_a.a + s->idx_n.a[rid], cnt, s->srt_a.a, s->srt_a.n, s->mini_cut, s->chain_cut, NULL);
// if(rid == 57) {
// fprintf(stderr, "-1-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
///remove candidate chains that have been calculated
if(!fi) backward_dedup_ol(rid, bl, &(b->sp), &b->olist);///it is ok
// if(rid == 57) {
// fprintf(stderr, "-2-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
filter_by_reliable_ovlp_mmhap_adv(rid, s->filter, &(b->sp), &b->olist, &(s->udb), s->sec_cutoff, 1, 1, s->ccov, &ol_h);
clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2);
if(!fi) ol_h = 0;
// if(rid == 57) {
// fprintf(stderr, "-3-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
ol_h = split_ug_lalign(ol_h, &b->olist, err_high, err_low,
&b->clist, &(s->udb), s->uopt, seq, len, &b->self_read, &b->ovlp_read,
&b->correct, &b->exz, aux_o, rid, s->k, s->chain_cut, NULL);
// if(rid == 57) {
// fprintf(stderr, "-4-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
aux_o = gen_aux_ovlp(&b->olist);///must be here
// if(rid == 57) {
// fprintf(stderr, "-5-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
ol_h = split_ug_lalign(ol_h, &b->olist, err_high, err_low,
&b->clist, &(s->udb), s->uopt, seq, len, &b->self_read, &b->ovlp_read,
&b->correct, &b->exz, aux_o, rid, s->k, s->chain_cut, NULL);
// if(rid == 57) {
// fprintf(stderr, "-6-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
if(fi) {///first round
pass_aln = test_het_aln_mmhap(rid, s->ccov, u_trans_a((*(s->filter)), rid), u_trans_n((*(s->filter)), rid), &b->olist, &(b->sp));
push_ul_ov_t(&(s->udb), u_trans_a((*(s->filter)), rid), u_trans_n((*(s->filter)), rid), rid, &(b->sp), &b->olist, len, pass_aln, err_high, bl);
// fprintf(stderr, "-1-[M::%s] utg%.6lu%c, rid::%lu, pass_aln::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, pass_aln);
} else {///second round
push_ul_ov_t(&(s->udb), NULL, 0, rid, &(b->sp), &b->olist, len, 0, err_high, bl);
remove_trans_ovlp_connect(s->udb.ug, rid, bl);
}
return pass_aln;
}
**/
static void worker_for_trans_chain_mmhap_adv(void *data, long i, int tid) // callback for kt_for()
{
ug_trans_t *s = (ug_trans_t*)data;
ha_ovec_buf_t *b = s->hab[tid]; kv_ul_ov_t *bl = &(s->ll[tid].tk);
uint32_t high_occ = asm_opt.polyploidy + 1; uint64_t cnt;
char *seq = s->ug->u.a[i].s; int64_t len = s->ug->u.a[i].len;
if((!s->is_ovlp) && (s->is_cnt)) s->idx_n.a[i] = 0;
if(s->ug->g->seq[i].del) return;
if(is_mmhom_node(s->ccov->cov.a+s->ccov->idx[i], &(s->ug->u.a[i]), s->ccov->rg, s->ccov->hom_min, 0.9)) return;
// asprintf(&as, "\n[M::%s] rid::%ld, len::%lu, name::%.*s\n", __func__, s->id+i, s->len[i], (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a);
// push_vlog(&(overall_zdbg->a[s->id+i]), as); free(as); as = NULL;
// if(rescue_adject_ovlp(s->ug->g, i, s->filter, &(s->ll[tid].lo))) return;
// gen_trans_chain_mmhap(s, i, b, bl, seq, len, 0.8, 0.8);
if(!s->is_ovlp) {
if(s->is_cnt) {
s->idx_n.a[i] = ug_map_lchain(b->abl, i, seq, len, s->w, s->k, &(s->udb), NULL, NULL, s->bw_thres, s->bw_thres_double,
s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, s->is_HPC, NULL, 0, NULL, 0, s->mini_cut, s->chain_cut, NULL);
} else {
cnt = ug_map_lchain(b->abl, i, seq, len, s->w, s->k, &(s->udb), NULL, NULL, s->bw_thres, s->bw_thres_double,
s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, s->is_HPC, s->idx_a.a + s->idx_n.a[i], 0, NULL, 0, s->mini_cut, s->chain_cut, NULL);
assert(cnt == ((s->idx_n.a[i+1]-s->idx_n.a[i])));
}
if(s->free_cnt[tid] >= FREE_BATCH) {
clear_count_buf(s, tid, 1); s->free_cnt[tid] = 0;
}
s->free_cnt[tid]++;
return;
}
// if(i == 58) {
// fprintf(stderr, "\n-1-[M::%s] utg%.6u%c, rid::%ld, is_ovlp::%d, is_cnt::%d, len::%ld, str::%u\n",
// __func__, (uint32_t)i+1, "lc"[s->ug->u.a[i].circ], i, s->is_ovlp, s->is_cnt, len, (uint32_t)(!!seq));
// }
if(!gen_trans_adaptive_mmhap_aln(s, i, b, bl, seq, len, s->filter, s->diff_ec_ul, s->diff_ec_ul_double, s->bw_thres, s->bw_thres_double)) {
gen_trans_adaptive_mmhap_aln(s, i, b, bl, seq, len, NULL, s->diff_ec_ul_double, s->diff_ec_ul_double, s->bw_thres_double, s->bw_thres_double);
}
if(s->free_cnt[tid] >= FREE_BATCH) {
clear_count_buf(s, tid, 0); s->free_cnt[tid] = 0;
}
s->free_cnt[tid]++;
}
int64_t retrieve_cigar_err_dir(bit_extz_t *ez, int64_t s, int64_t e, int64_t *xk, int64_t *ck, int64_t is_back)
{ ///[ez->ts, ez->te]/[ez->qs, ez->qe]/[s, e)
@@ -20178,6 +20374,145 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res)
fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time);
}
void clean_trans_base_count_mmhap_comp_rmap(ug_trans_t *p, kv_u_trans_t *res)
{
uint64_t i, k, l, occ, idx_n; ha_mzl_t *tz; u_trans_t *z;
kv_ul_ov_t *bl; double ww; ha_mzl_t *idx;
///make results consistent
kv_resize(ha_mzl_t, p->srt_a, p->srt_a.n+p->ug->u.n);
idx = p->srt_a.a + p->srt_a.n; idx_n = p->ug->u.n;
for (i = 0; i < idx_n; i++) {
tz = &(idx[i]);
tz->x = (uint64_t)-1; tz->rev = 0;
tz->pos = tz->rid = tz->span = 0;
}
for (i = 0, occ = res->n; (int64_t)i < p->n_thread; i++) {
bl = &(p->ll[i].tk);
if(!(bl->n)) continue;
for (k = 1, l = 0; k <= bl->n; k++) {
if(k == bl->n || bl->a[k].qn != bl->a[l].qn) {
if(k > l) {
tz = &(idx[bl->a[l].qn]);
tz->x = bl->a[l].qn; tz->x <<= 32; tz->x |= i;
tz->rid = l>>32; tz->pos = (uint32_t)l; tz->rev = 1;
occ += (k - l);
}
l = k;
}
}
}
kv_resize(u_trans_t, *res, occ);
for (i = 0; i < idx_n; i++) {
tz = &(idx[i]);
if(!(tz->rev)) continue;
bl = &(p->ll[(uint32_t)(tz->x)].tk);
k = tz->rid; k <<= 32; k += tz->pos;
assert(bl->a[k].qn == (tz->x>>32));
for (; (k < bl->n) && (bl->a[k].qn == (tz->x>>32)); k++) {
if(bl->a[k].qn == bl->a[k].tn) continue;
ww = cal_trans_ov_w(&(bl->a[k]));
if(ww <= 0) continue;
kv_pushp(u_trans_t, *res, &z);
z->f = RC_3; z->rev = bl->a[k].rev; z->del = 0;
z->qn = bl->a[k].qn; z->qs = bl->a[k].qs; z->qe = bl->a[k].qe;
z->tn = bl->a[k].tn; z->ts = bl->a[k].ts; z->te = bl->a[k].te;
z->nw = ww;
}
}
destory_ug_rid_cov_t(p->ccov); free(p->ccov);
p->ccov = gen_ug_rid_cov_t(p->ug, p->rg, p->uopt->sources);
clean_u_trans_t_idx_filter_mmhap_adv(res, p->ug, p->rg, p->uopt->sources, p->ccov);
gen_ug_rid_cov_t_by_ovlp(res, p->ccov);
}
void gen_trans_base_count_mmhap_comp_rmap(ug_trans_t *p, kv_u_trans_t *res)
{
double index_time = yak_realtime();
uint64_t i, k, l, occ, m, cc;
p->ccov = gen_ug_rid_cov_t(p->ug, p->rg, p->uopt->sources);
clean_u_trans_t_idx_adv(res, p->ug, p->rg); p->filter = res;
p->is_cnt = 1; p->is_ovlp = 0;
memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread);
kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n);
for (i = l = 0; i < p->ug->u.n; i++) {
occ = p->idx_n.a[i]; p->idx_n.a[i] = l; l += occ;
}
p->idx_n.a[i] = l;
p->idx_a.n = p->idx_a.m = l; MALLOC(p->idx_a.a, p->idx_a.n);
p->is_cnt = 0; p->is_ovlp = 0;
memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread);
kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n);
p->srt_a.n = p->srt_a.m = p->idx_a.n; MALLOC(p->srt_a.a, p->srt_a.n);
for (i = 0; i < p->srt_a.n; i++) {
p->srt_a.a[i] = p->idx_a.a[i];
p->srt_a.a[i].pos = (uint32_t)i;
p->srt_a.a[i].rid = i>>32;
}
radix_sort_ha_mzl_t_srt(p->srt_a.a, p->srt_a.a + p->srt_a.n);
kvec_t(uint64_t) cut; kv_init(cut);
for (k = 1, l = 0; k <= p->srt_a.n; k++) {
if(k == p->srt_a.n || p->srt_a.a[l].x != p->srt_a.a[k].x) {
for (i = l; i < k; i++) {
m = p->srt_a.a[i].rid; m <<= 32; m |= p->srt_a.a[i].pos;
assert(p->srt_a.a[i].x == p->idx_a.a[m].x);
p->srt_a.a[i] = p->idx_a.a[m]; p->idx_a.a[m].x = i;
}
kv_push(uint64_t, cut, (k - l));
l = k;
}
}
if(cut.n > 0) {
radix_sort_gfa64(cut.a, cut.a + cut.n);
m = cut.n * 0.0002; cc = cut.a[cut.n-1] + 1;
if(m > 0 && m <= cut.n) cc = cut.a[cut.n-m] + 1;
if(cc < (uint64_t)p->mini_cut) p->mini_cut = cc;
}
kv_destroy(cut);
p->is_cnt = 0; p->is_ovlp = 1;
memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread);
kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n);
clean_trans_base_count_mmhap_comp_rmap(p, res);
for (i = 0; (int64_t)i < p->n_thread; i++) p->ll[i].tk.n = 0;
p->is_cnt = 0; p->is_ovlp = 1;
memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread);
kt_for(p->n_thread, worker_for_trans_chain_mmhap_adv, p, p->ug->u.n);
for (i = 0; (int64_t)i < p->n_thread; i++) {
ha_ovec_destroy(p->hab[i]);
free(p->ll[i].lo.a); free(p->ll[i].srt.a.a); free(p->ll[i].tc.a);
}
free(p->idx_a.a); free(p->idx_n.a); free(p->hab); free(p->free_cnt);
destory_ug_rid_cov_t(p->ccov); free(p->ccov);
for (i = 0; (int64_t)i < p->n_thread; i++) free(p->ll[i].tk.a);
free(p->srt_a.a); free(p->ll);
fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time);
}
void gen_trans_base_count_mmhap_comp(ug_trans_t *p, kv_u_trans_t *res)
{