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

backup r321

Browse Source
This commit is contained in:
chhylp123
2021-04-01 21:02:57 -04:00
parent 0605aa1961
commit 46e899bbae
11 changed files with 1494 additions and 401 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+16 -1
View File
@@ -32,6 +32,8 @@ static ko_longopt_t long_options[] = {
{ "b-cov", ko_required_argument, 317 },
{ "h-cov", ko_required_argument, 318 },
{ "m-rate", ko_required_argument, 319 },
{ "b-partition", ko_no_argument, 320 },
{ "t-occ", ko_required_argument, 321 },
{ 0, 0, 0 }
};
@@ -87,6 +89,9 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " -d INT upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt);
fprintf(stderr, " -3 FILE list of hap1/paternal read names []\n");
fprintf(stderr, " -4 FILE list of hap2/maternal read names []\n");
fprintf(stderr, " --t-occ INT\n");
fprintf(stderr, " force remove unitigs with >INT unexpected haplotype-specific reads;\n");
fprintf(stderr, " ignore graph topology; [%d]\n", asm_opt->trio_flag_occ_thres);
fprintf(stderr, " Purge-dups:\n");
fprintf(stderr, " -l INT purge level. 0: no purging; 1: light; 2/3: aggressive [0 for trio; 2 for unzip]\n");
@@ -96,6 +101,7 @@ void Print_H(hifiasm_opt_t* asm_opt)
asm_opt->purge_overlap_len);
fprintf(stderr, " --purge-cov INT\n");
fprintf(stderr, " coverage upper bound of Purge-dups [auto]\n");
///fprintf(stderr, " --high-het enable this mode for high heterozygosity sample [experimental, not stable]\n");
fprintf(stderr, " Hi-C-partition [experimental, not stable]:\n");
@@ -157,6 +163,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->recover_atg_cov_min = -1024;
asm_opt->recover_atg_cov_max = INT_MAX;
asm_opt->hom_global_coverage = -1;
asm_opt->hom_global_coverage_set = 0;
asm_opt->bed_inconsist_rate = 70;
asm_opt->hic_inconsist_rate = 30;
///asm_opt->bub_mer_length = 3;
@@ -165,6 +172,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->b_high_cov = -1;
asm_opt->m_rate = 0.75;
asm_opt->hap_occ = 1;
asm_opt->polyploidy = 2;
asm_opt->trio_flag_occ_thres = 60;
}
void destory_enzyme(enzyme* f)
@@ -616,7 +625,11 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 306) asm_opt->max_ov_diff_final = atof(opt.arg);
else if (c == 307) asm_opt->extract_list = opt.arg;
else if (c == 308) asm_opt->extract_iter = atoi(opt.arg);
else if (c == 309) asm_opt->hom_global_coverage = atoi(opt.arg);
else if (c == 309)
{
asm_opt->hom_global_coverage = atoi(opt.arg);
asm_opt->hom_global_coverage_set = 1;
}
else if (c == 310)
{
char* s = NULL;
@@ -636,6 +649,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 317) asm_opt->b_low_cov = atoi(opt.arg);
else if (c == 318) asm_opt->b_high_cov = atoi(opt.arg);
else if (c == 319) asm_opt->m_rate = atof(opt.arg);
else if (c == 320) asm_opt->flag |= HA_F_PARTITION;
else if (c == 321) asm_opt->trio_flag_occ_thres = atoi(opt.arg);
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);
CommandLines.h
+5 -1
View File
@@ -3,7 +3,7 @@
#include <pthread.h>
#define HA_VERSION "0.14.2-r317"
#define HA_VERSION "0.14.2-r320"
#define VERBOSE 0
@@ -18,6 +18,7 @@
#define HA_F_BAN_POST_JOIN 0x100
#define HA_F_BAN_ASSEMBLY 0x200
#define HA_F_HIGH_HET 0x400
#define HA_F_PARTITION 0x800
#define HA_MIN_OV_DIFF 0.02 // min sequence divergence in an overlap
@@ -74,6 +75,7 @@ typedef struct {
int recover_atg_cov_min;
int recover_atg_cov_max;
int hom_global_coverage;
int hom_global_coverage_set;
int bed_inconsist_rate;
int hic_inconsist_rate;
@@ -93,6 +95,8 @@ typedef struct {
long long mem_buf;
long long coverage;
int hap_occ;
int polyploidy;
int trio_flag_occ_thres;
} hifiasm_opt_t;
Makefile
+2 -1
View File
@@ -6,7 +6,7 @@ CPPFLAGS=
INCLUDES=
OBJS= CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \
htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o
htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o partig.o
EXE= hifiasm
LIBS= -lz -lpthread -lm
@@ -72,3 +72,4 @@ main.o: Levenshtein_distance.h htab.h
sketch.o: kvec.h htab.h Process_Read.h Overlaps.h kdq.h CommandLines.h
sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
hic.o: hic.h
partig.o: partig.h
Overlaps.cpp
+236 -89
View File
@@ -45,7 +45,8 @@ KSORT_INIT_GENERIC(uint32_t)
///this value has been updated at the first line of build_string_graph_without_clean
long long min_thres;
uint32_t print_untig_by_read(ma_ug_t *g, const char* name, uint32_t in, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, const char* info);
void init_bub_label_t(bub_label_t* x, uint32_t n_thres, uint32_t n_reads)
{
@@ -11711,61 +11712,6 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp)
kv_destroy(new_rtg_edges.a);
}
void classify_untigs_debug(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex,
kvec_asg_arc_t_warp* new_rtg_edges, int max_hang, int min_ovlp, trans_chain* t_ch)
{
uint64_t i, dip_thre_max, dip_thres, n_utg;
uint8_t* primary_flag = (uint8_t*)calloc(sg->n_seq, sizeof(uint8_t));
hap_cov_t *cov = init_hap_cov_t(ug, sg, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 0);
int is_set = ((asm_opt.hom_global_coverage != -1)? 1 : 0);
if(is_set == 0)
{
purge_dups(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges,
asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, 0, 0, 1, cov, 0);
dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE);
asm_opt.hom_global_coverage = -1;
}
else
{
dip_thre_max = asm_opt.hom_global_coverage;
}
dip_thre_max *= 0.75;
for (i = 0; i < ug->g->n_seq; i++)
{
dip_thres = dip_thre_max;
if(ug->u.a[i].n <= dip_thre_max) dip_thres = dip_thre_max * 0.9;
if(get_ug_coverage(&ug->u.a[i], sg, coverage_cut, sources, ruIndex, primary_flag)<dip_thres)
{
ug->g->seq[i].c = 1;
if(get_ug_coverage_aggressive(ug, i, sg, coverage_cut, sources, ruIndex, primary_flag, &n_utg)<dip_thres)
{
if(n_utg > dip_thres*2) ug->g->seq[i].c = 2;
}
}
else
{
ug->g->seq[i].c = 0;
}
uint32_t flag = get_unitig_het_arb(&ug->u.a[i], t_ch->is_r_het, 2, 1, 0);
if(flag != ug->g->seq[i].c)
{
if(flag == 0 || ug->g->seq[i].c == 0) fprintf(stderr, "Attention\n");
fprintf(stderr, "utg%.6lul, flag: %u, ug->g->seq[i].c: %u\n", i+1, flag, ug->g->seq[i].c);
}
}
free(primary_flag);
destory_hap_cov_t(&cov);
fprintf(stderr, "[M::%s] diploid coverage threshold: %lu\n", __func__, dip_thre_max);
exit(0);
}
void set_ug_coverage_aggressive(ma_ug_t *ug, uint32_t uID, asg_t* read_g,
const ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag,
trans_chain* t_ch, long long het_cov_thres)
@@ -12203,6 +12149,180 @@ bub_label_t* b_mask_t)
0.05, 0.9, max_hang, min_ovlp, 0, b_mask_t);
}
void set_trio_flag_by_cov(ma_ug_t *ug, hap_cov_t *cov)
{
kvec_t(uint64_t) idx; kv_init(idx);
uint32_t i, k, j, qn, tn, s[2], flag;
ma_utg_t *u = NULL, *w = NULL;
hc_links link;
init_hc_links(&link, ug->g->n_seq, cov->t_ch);
for (i = 0, idx.n = 0; i < link.a.n; i++)
{
qn = i;
u = &(ug->u.a[qn]);
for (k = 0; k < u->n; k++)
{
if((R_INF.trio_flag[u->a[k]>>33]&SET_TRIO)==0) break;
}
if(k >= u->n) continue; ///whole unitig is primary
for (k = 0, s[0] = s[1] = 0; k < link.a.a[qn].f.n; k++)
{
tn = link.a.a[qn].f.a[k].uID;
w = &(ug->u.a[tn]);
for (j = 0; j < w->n; j++)
{
if((R_INF.trio_flag[w->a[j]>>33]&FATHER)||(R_INF.trio_flag[w->a[j]>>33]&MOTHER))
{
s[0]++;
}
if(R_INF.trio_flag[w->a[j]>>33]&SET_TRIO)
{
s[1]++;
}
}
}
if(s[1] > 0)
{
kv_push(uint64_t, idx, (uint64_t)((uint32_t)-1 - s[0]) << 32 | (qn));
}
}
for (i = 0; i < idx.n; i++)
{
qn = (uint32_t)idx.a[i];
u = &(ug->u.a[qn]);
for (k = 0, s[0] = s[1] = 0; k < link.a.a[qn].f.n; k++)
{
tn = link.a.a[qn].f.a[k].uID;
w = &(ug->u.a[tn]);
for (j = 0; j < w->n; j++)
{
if(R_INF.trio_flag[w->a[j]>>33]&FATHER)
{
s[0]++;
}
if(R_INF.trio_flag[w->a[j]>>33]&MOTHER)
{
s[1]++;
}
}
}
if(s[0] >= s[1])
{
flag = MOTHER;
}
else
{
flag = FATHER;
}
for (k = 0; k < u->n; k++)
{
if(R_INF.trio_flag[u->a[k]>>33]&SET_TRIO) continue;
if(cov->t_ch->is_r_het[u->a[k]>>33] == N_HET) continue;
R_INF.trio_flag[u->a[k]>>33] |= flag;
}
}
for (i = 0, idx.n = 0; i < link.a.n; i++)
{
qn = i;
u = &(ug->u.a[qn]);
for (k = 0; k < u->n; k++)
{
if(R_INF.trio_flag[u->a[k]>>33]&FATHER)
{
R_INF.trio_flag[u->a[k]>>33] = FATHER;
}
else if(R_INF.trio_flag[u->a[k]>>33]&MOTHER)
{
R_INF.trio_flag[u->a[k]>>33] = MOTHER;
}
else
{
R_INF.trio_flag[u->a[k]>>33] = AMBIGU;
}
}
}
kv_destroy(idx);
destory_hc_links(&link);
}
void print_r_het(hap_cov_t *cov, uint8_t* trio_flag, const char* cmd)
{
if(cov && cov->t_ch)
{
fprintf(stderr, "\n+%s-is_r_het[1369536]=%u\n", cmd, cov->t_ch->is_r_het[1369536]);
fprintf(stderr, "+%s-is_r_het[5097804]=%u\n", cmd, cov->t_ch->is_r_het[5097804]);
fprintf(stderr, "+%s-is_r_het[603738]=%u\n", cmd, cov->t_ch->is_r_het[603738]);
}
if(trio_flag)
{
fprintf(stderr, "+%s-trio_flag[1369536]=%u\n", cmd, trio_flag[1369536]);
fprintf(stderr, "+%s-trio_flag[5097804]=%u\n", cmd, trio_flag[5097804]);
fprintf(stderr, "+%s-trio_flag[603738]=%u\n", cmd, trio_flag[603738]);
}
}
void output_bp_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
long long tipsLen, float tip_drop_ratio, long long stops_threshold,
R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
bub_label_t* b_mask_t)
{
hic_clean(sg);
kvec_asg_arc_t_warp new_rtg_edges;
kv_init(new_rtg_edges.a);
ma_ug_t *ug = NULL;
ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
hap_cov_t *cov = NULL;
asg_t *copy_sg = copy_read_graph(sg);
ma_ug_t *copy_ug = copy_untig_graph(ug);
///asm_opt.purge_overlap_len = asm_opt.purge_overlap_len_hic;
///asm_opt.purge_simi_thres = asm_opt.purge_simi_rate_hic;
adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut,
tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio,
max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1);
// print_untig_by_read(copy_ug, "m64011_190830_220126/175638789/ccs", 1369536, NULL, NULL, "sa");
// print_untig_by_read(copy_ug, "m64012_190921_234837/21039588/ccs", 5097804, NULL, NULL, "sa");
// print_untig_by_read(copy_ug, "m64011_190830_220126/88867583/ccs", 603738, NULL, NULL, "sa");
ma_ug_destroy(copy_ug);
asg_destroy(copy_sg);
// print_untig_by_read(copy_ug, "m64011_190830_220126/175638789/ccs", 1369536, NULL, NULL, "sb");
// print_untig_by_read(copy_ug, "m64012_190921_234837/21039588/ccs", 5097804, NULL, NULL, "sb");
// print_untig_by_read(copy_ug, "m64011_190830_220126/88867583/ccs", 603738, NULL, NULL, "sb");
// print_r_het(cov, R_INF.trio_flag, "out-0");
set_trio_flag_by_cov(ug, cov);
// print_r_het(cov, R_INF.trio_flag, "out-1");
destory_hap_cov_t(&cov);
ma_ug_destroy(ug);
kv_destroy(new_rtg_edges.a);
output_unitig_graph(sg, coverage_cut, output_file_name, sources, ruIndex, max_hang, min_ovlp);
output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex,
0.05, 0.9, max_hang, min_ovlp, 0, b_mask_t);
output_trio_unitig_graph(sg, coverage_cut, output_file_name, MOTHER, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex,
0.05, 0.9, max_hang, min_ovlp, 0, b_mask_t);
}
ma_ug_t* merge_utg(ma_ug_t **dest, ma_ug_t **src)
{
asg_t *g_d = (*dest)->g, *g_s = (*src)->g;
@@ -13039,7 +13159,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre
}
}
if(cov && operation != CUT) collect_trans_cov(__func__, &b_0, &b_1, ug, read_sg, cov);
///note: we need to remove b_0, insetad of b_1 here
if(cov && operation != CUT) collect_trans_cov(__func__, &b_1, &b_0, ug, read_sg, cov);
break;
}
@@ -13388,8 +13509,6 @@ R_to_U* ruIndex, uint32_t positive_flag, float drop_rate)
}
}
del_node++;
}
@@ -13517,7 +13636,8 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_
asg_seq_drop(g, b.b.a[k]>>1);
}
if(cov) collect_trans_cov(__func__, &b_0, &b_1, ug, read_sg, cov);
///note: we need to remove b_0, insetad of b_1 here
if(cov) collect_trans_cov(__func__, &b_1, &b_0, ug, read_sg, cov);
///lable the primary one
b_0.b.n = 0;
@@ -13960,7 +14080,7 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, trio_flag, drop_ratio);
drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, DOUBLE_CHECK_THRES);
all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, DOUBLE_CHECK_THRES, asm_opt.trio_flag_occ_thres);
if(is_first)
{
is_first = 0;
@@ -13986,6 +14106,8 @@ void print_graph_statistic(asg_t *g, const char* cmd)
fprintf(stderr, "%s->n_node: %lu, n_arc: %lu, size: %lu\n", cmd, n_node, n_arc, size);
}
void clean_primary_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, long long tipsLen, float tip_drop_ratio,
long long stops_threshold, R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density,
@@ -13997,18 +14119,16 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov)
int round = T_ROUND;
///uint32_t is_first = 1;
redo:
redo:
///print_graph_statistic(g, "beg");
///print_debug_gfa(read_g, ug, coverage_cut, "debug_trans_ovlp_hg002", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, /**is_first?0:1**/1);
///untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, DROP, cov);
///is_first = 0;
if(just_bubble_pop == 0)
{
cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2);
}
long long pre_cons = get_graph_statistic(g);
long long cur_cons = 0;
while(pre_cons != cur_cons)
@@ -14038,7 +14158,6 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov)
cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex,
2);
}
resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, (uint32_t)-1, drop_ratio);
drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, reverse_sources, 0, 1);
@@ -14075,8 +14194,6 @@ void set_drop_trio_flag(ma_ug_t *ug)
}
}
void update_unitig_graph(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex,
uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate)
@@ -14109,6 +14226,8 @@ uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate)
if(R_INF.trio_flag[rId] == flag) flag_occ++;
if(R_INF.trio_flag[rId] != flag) non_flag_occ++;
}
///if(is_final_check && v == 0) fprintf(stderr, "flag: %u, flag_occ: %u, non_flag_occ: %u\n", flag, flag_occ, non_flag_occ);
if(is_final_check == 0 && hap_label_occ == u->n) continue;
///if(is_double_check && non_flag_occ < u->n*DOUBLE_CHECK_THRES) continue;
@@ -14169,7 +14288,6 @@ uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate)
}
void get_candidate_uids(asg_t* nsg, ma_utg_t* nsu, kvec_t_u64_warp* u_vecs,
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex)
{
@@ -14290,13 +14408,13 @@ uint32_t* non_require, uint32_t* ambigious)
///note: to use this function, don't renew unitig graph!!!!!!!!!
void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate,
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate)
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate, int non_tig_occ)
{
kvec_t_u64_warp u_vecs;
kv_init(u_vecs.a);
uint32_t n_vtx, v, k, is_Unitig, uId, rId, convex, flag_occ, non_flag_occ, ambigious, /**is_ambigious,**/ flag, n_reduce = 1;
uint32_t n_vtx, v, k, is_Unitig, is_tig, uId, rId, convex, flag_occ, non_flag_occ, ambigious, /**is_ambigious,**/ flag, n_reduce = 1;
asg_t* nsg = NULL;
ma_utg_t* nsu = NULL;
nsg = ug->g;
@@ -14332,11 +14450,16 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate)
if(nsu->m == 0) continue;
if(nsg->seq[uId].del) continue;
if(nsg->seq[uId].c == ALTER_LABLE) continue;
is_tig = 1;
if(get_real_length(nsg, v^1, NULL) == 1)
{
get_real_length(nsg, v^1, &convex);
if(get_real_length(nsg, convex^1, NULL) == 1) continue;
if(get_real_length(nsg, convex^1, NULL) == 1)
{
is_tig = 0;
///continue;
}
}
get_unitig_trio_flag(nsu, flag, &flag_occ, &non_flag_occ, &ambigious);
@@ -14349,6 +14472,7 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate)
else**/
{
///we may need it or not
if(is_tig == 0 && (int)flag_occ <= non_tig_occ) continue;
if(flag_occ <= ((non_flag_occ+flag_occ)*drop_rate)) continue;
if((flag_occ+non_flag_occ) == 0) continue;
if(nsu->n >= 100)
@@ -14698,7 +14822,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
uint32_t v, n_vtx = nsg->n_seq;
hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources,
coverage_cut, max_hang, min_ovlp, asm_opt.purge_level_trio>0?1:0);
if(asm_opt.purge_level_primary > 0)
if(cov->t_ch)
{
set_r_het_flag(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex,
new_rtg_edges, max_hang, min_ovlp, cov->t_ch);
@@ -14742,6 +14866,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
clean_trio_untig_graph(*ug, read_g, coverage_cut, sources, reverse_sources, tipsLen,
tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio, flag, drop_rate, cov);
///delete_useless_nodes(ug);
delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex);
@@ -14750,6 +14875,8 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0,
DOUBLE_CHECK_THRES, flag, drop_rate);
///if(flag == MOTHER) print_debug_gfa(read_g, *ug, coverage_cut, "debug_trio_1", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
renew_utg(ug, read_g, new_rtg_edges);
@@ -14767,7 +14894,8 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
}
///if(flag == MOTHER) print_untig_by_read(*ug, "m64043_200627_000137/124716590/ccs", 2789716, NULL, NULL, "beg");
update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 1,
FINAL_DOUBLE_CHECK_THRES, flag, drop_rate);
@@ -22502,8 +22630,8 @@ uint32_t collect_p_trans)
uint32_t v, n_vtx = nsg->n_seq, k, rId, just_contain;
ma_utg_t* u = NULL;
hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources,
coverage_cut, max_hang, min_ovlp, asm_opt.purge_level_primary>0?1:0);
if(asm_opt.purge_level_primary > 0)
coverage_cut, max_hang, min_ovlp, (asm_opt.purge_level_primary>0||i_cov)?1:0);
if(cov->t_ch)
{
set_r_het_flag(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex,
new_rtg_edges, max_hang, min_ovlp, cov->t_ch);
@@ -22519,13 +22647,13 @@ uint32_t collect_p_trans)
nsg->seq[v].c = PRIMARY_LABLE;
EvaluateLen((*ug)->u, v) = (*ug)->u.a[v].n;
}
clean_primary_untig_graph(*ug, read_g, sources, reverse_sources, coverage_cut, tipsLen, tip_drop_ratio,
stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio, cov);
delete_useless_nodes(ug);
renew_utg(ug, read_g, new_rtg_edges);
if(i_cov) goto skip_purge;
///classify_untigs_debug(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, max_hang, min_ovlp, cov->t_ch);
if(asm_opt.purge_level_primary > 0)
{
@@ -22620,16 +22748,13 @@ uint32_t collect_p_trans)
recover_utg_by_coverage(ug, read_g, coverage_cut, sources, ruIndex, cov->t_ch);
if(i_cov)
{
if(asm_opt.purge_level_primary > 0)
if(cov->t_ch && collect_p_trans)
{
if(collect_p_trans)
{
// purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges,
// asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio,
// 0, 0, cov, collect_p_trans);
}
purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges,
asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio,
0, 0, cov, collect_p_trans);
}
(*i_cov) = cov;
}
else
@@ -26624,7 +26749,6 @@ bub_label_t* b_mask_t)
asg_destroy(copy_sg); copy_sg = NULL;
uint32_t beg_idx, occ;
bubble_type bub;
memset(&bub, 0, sizeof(bubble_type));
@@ -26634,15 +26758,12 @@ bub_label_t* b_mask_t)
rescue_bubbles_by_contained_reads(ug, sg, sources, coverage_cut, ruIndex, max_hang, min_ovlp, chainLenThres, beg_idx, occ, &bub, b_mask_t);
///output_unitig_graph(sg, coverage_cut, (char*)"debug_1.rescue", sources, ruIndex, max_hang, min_ovlp);
ma_ug_destroy(ug); ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
reset_bub(&bub, ug, cov->t_ch, &new_rtg_edges);
beg_idx = bub.f_bub; occ = bub.b_bub + bub.b_end_bub + bub.tangle_bub;
rescue_bubbles_by_missing_ovlp(ug, sg, sources, coverage_cut, ruIndex, max_hang, min_ovlp, chainLenThres, beg_idx, occ, &bub, b_mask_t);
///output_unitig_graph(sg, coverage_cut, (char*)"debug_2.hic", sources, ruIndex, max_hang, min_ovlp);
ma_ug_destroy(ug); ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
reset_bub(&bub, ug, cov->t_ch, &new_rtg_edges);
beg_idx = bub.f_bub; occ = bub.b_bub + bub.b_end_bub + bub.tangle_bub;
@@ -26655,7 +26776,6 @@ bub_label_t* b_mask_t)
reset_bub(&bub, ug, cov->t_ch, &new_rtg_edges);
rescue_missing_hap_ovlp(ug, sg, sources, coverage_cut, max_hang, min_ovlp, &bub, gap_fuzz);
}
destory_bubbles(&bub);
@@ -28490,6 +28610,22 @@ int max_hang, int min_ovlp)
}
}
void set_hom_global_coverage(hifiasm_opt_t *opt, asg_t *sg, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, int max_hang, int min_ovlp)
{
if(opt->hom_global_coverage_set == 0)
{
ma_ug_t *ug = NULL;
ug = ma_ug_gen(sg);
hap_cov_t *cov = init_hap_cov_t(ug, sg, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 0);
purge_dups(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, NULL,
opt->purge_simi_thres, opt->purge_overlap_len, max_hang, min_ovlp, 0, 0, 1, cov, 0);
destory_hap_cov_t(&cov);
ma_ug_destroy(ug);
}
}
void clean_graph(
int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
@@ -28694,6 +28830,9 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
/*******************************for debug***************************************/
}
// set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex,
// max_hang_length, mini_overlap_length);
if (ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt))
{
@@ -28725,6 +28864,14 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, &b_mask_t);
free(buf);
}
else if(asm_opt.flag & HA_F_PARTITION)
{
char *buf = (char*)calloc(strlen(output_file_name) + 25, 1);
sprintf(buf, "%s.bp", output_file_name);
output_bp_graph(sg, coverage_cut, buf, sources, reverse_sources, (asm_opt.max_short_tip*2),
0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, &b_mask_t);
free(buf);
}
else
{
output_unitig_graph(sg, coverage_cut, output_file_name, sources, ruIndex, max_hang_length, mini_overlap_length);
Overlaps.h
+2 -1
View File
@@ -1013,7 +1013,7 @@ uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp*
void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t is_bubble_check, uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t);
void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate);
ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate, int non_tig_occ);
void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges, bub_label_t* b_mask_t);
@@ -1056,6 +1056,7 @@ typedef struct{
#define N_HET 0
#define C_HET 1
#define P_HET 2
#define S_HET 4
typedef struct{
kvec_t(uint32_t) uIDs;
Process_Read.h
+1
View File
@@ -101,6 +101,7 @@ typedef struct
#define MIX_TRIO 3
#define NON_TRIO 4
#define DROP 5
#define SET_TRIO 8
typedef struct
{
Purge_Dups.cpp
+196 -294
View File
@@ -8,6 +8,7 @@
#include "kthread.h"
#include "kdq.h"
#include "hic.h"
#include "partig.h"
KDQ_INIT(uint64_t)
@@ -49,33 +50,6 @@ typedef struct {
#define MIXED 3
#define FLIP 4
typedef struct {
uint8_t rev;
uint8_t type;
uint8_t status;
uint32_t x_beg_pos;
uint32_t x_end_pos;
uint32_t y_beg_pos;
uint32_t y_end_pos;
uint32_t x_beg_id;
uint32_t x_end_id;
uint32_t y_beg_id;
uint32_t y_end_id;
uint32_t xUid;
uint32_t yUid;
uint32_t weight;
long long score;
}hap_overlaps;
typedef struct {
kvec_t(hap_overlaps) a;
}kvec_hap_overlaps;
typedef struct {
kvec_hap_overlaps* x;
uint32_t num;
}hap_overlaps_list;
typedef struct {
uint64_t* vote_counting;
@@ -2377,7 +2351,7 @@ ma_hit_t_alloc* reverse_sources, long long xBegPos, long long xEndPos)
// fprintf(stderr, "tailIndex->a.n: %u, xReads->n: %u, total_match: %lld, hap_match: %lld, inp_match: %lld\n",
// tailIndex->a.n, xReads->n, (xEndPos - xBegPos + 1), hap_match, inp_match);
return (inp_match*3) - ((hap_match-inp_match)*1);
return ((double)(inp_match)*1) - ((double)(hap_match-inp_match)*0.334);
}
@@ -3139,7 +3113,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid)
}
int inline get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn)
int get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn)
{
uint32_t i;
for (i = 0; i < x->a.n; i++)
@@ -3324,6 +3298,10 @@ double filter_rate)
{
homLen += ovlp;
}
else if(asm_opt.polyploidy <= 2 && (a[k].h_status&S_HET))
{
homLen += ovlp;
}
else
{
hetLen += ovlp;
@@ -3354,10 +3332,14 @@ double filter_rate)
{
homLen += ovlp;
}
else if(asm_opt.polyploidy <= 2 && (a[k].h_status&S_HET))
{
homLen += ovlp;
}
else
{
hetLen += ovlp;
}
}
}
if(hetLen <= ((hetLen + homLen) * filter_rate))
@@ -4774,242 +4756,6 @@ void remove_contained_haplotig(hap_overlaps_list* all_ovlp, ma_ug_t *ug, asg_t*
// }
}
#define HAP1_LABLE 7
#define HAP2_LABLE 8
void init_contig_phase(hap_overlaps_list* all_ovlp, asg_t *purge_g, float drop_ratio)
{
uint32_t v, i, n_vtx = purge_g->n_seq * 2, beg, end, uId, yId, need_update;
long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen, hap1_weight, hap2_weight;
buf_t b_0;
memset(&b_0, 0, sizeof(buf_t));
clean_purge_graph(purge_g, drop_ratio, 1);
for (v = 0; v < n_vtx; ++v)
{
if(purge_g->seq[v>>1].c == ALTER_LABLE || purge_g->seq[v>>1].del) continue;
if(get_real_length(purge_g, v, NULL) != 1) continue;
if(get_real_length(purge_g, v^1, NULL) != 0) continue;
beg = v;
b_0.b.n = 0;
if(get_unitig(purge_g, NULL, beg, &end, &nodeLen, &baseLen, &max_stop_nodeLen,
&max_stop_baseLen, 1, &b_0) == LOOP)
{
continue;
}
for (i = 0; i < b_0.b.n; i++)
{
uId = b_0.b.a[i];
if((i>>1) == 0) purge_g->seq[uId].c = HAP1_LABLE;
else purge_g->seq[uId].c = HAP2_LABLE;
}
}
need_update = 1;
while (need_update)
{
need_update = 0;
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
if(purge_g->seq[uId].c == HAP1_LABLE || purge_g->seq[uId].c == HAP2_LABLE)
{
continue;
}
if(all_ovlp->x[uId].a.n == 0)
{
purge_g->seq[uId].c = HAP1_LABLE;
continue;
}
hap1_weight = hap2_weight = 0;
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
yId = all_ovlp->x[uId].a.a[i].yUid;
if(purge_g->seq[yId].c == HAP1_LABLE) hap1_weight += all_ovlp->x[uId].a.a[i].weight;
if(purge_g->seq[yId].c == HAP2_LABLE) hap2_weight += all_ovlp->x[uId].a.a[i].weight;
}
if(hap1_weight == 0 && hap2_weight == 0)
{
need_update = 1;
continue;
}
if(hap1_weight >= hap2_weight)
{
purge_g->seq[uId].c = HAP1_LABLE;
}
else
{
purge_g->seq[uId].c = HAP2_LABLE;
}
}
}
free(b_0.b.a);
}
void partition_contigs(hap_overlaps_list* all_ovlp, ma_ug_t *ug, asg_t *purge_g, hap_cov_t *cov, double keep_rate,
int max_hang, int min_ovlp, float drop_ratio)
{
int r, index;
long long max_score;
uint32_t v, i, uId, xUid, is_contain = 0, m;
hap_overlaps *p = NULL;
asg_arc_t t, *p_t;
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(all_ovlp->x[uId].a.a[i].status == DELETE) continue;
all_ovlp->x[uId].a.a[i].status = MIXED;
}
}
for (v = 0; v < all_ovlp->num; v++)
{
uId = v; p = NULL; is_contain = 0;
if(all_ovlp->x[uId].a.n == 0) continue;
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(p == NULL || p->score < all_ovlp->x[uId].a.a[i].score)
{
p = &(all_ovlp->x[uId].a.a[i]);
}
}
max_score = p->score;
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(all_ovlp->x[uId].a.a[i].type == YCX)
{
if(!filter_secondary_chain(p->score, all_ovlp->x[uId].a.a[i].score, MAX(0.95, keep_rate)))
{
continue;
}
xUid = all_ovlp->x[uId].a.a[i].xUid;
purge_g->seq[xUid].c = ALTER_LABLE;
purge_g->seq[xUid].del = 1;
///collect_trans_purge_cov(cov, ug, &(all_ovlp->x[uId].a.a[i]), 0);
if(is_contain == 0 || all_ovlp->x[uId].a.a[i].score > max_score)
{
max_score = all_ovlp->x[uId].a.a[i].score;
}
is_contain = 1;
}
}
if(is_contain)
{
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(filter_secondary_chain(max_score, all_ovlp->x[uId].a.a[i].score, keep_rate))
{
all_ovlp->x[uId].a.a[i].status = FLIP;
index = get_specific_hap_overlap(&(all_ovlp->x[all_ovlp->x[uId].a.a[i].yUid]),
all_ovlp->x[uId].a.a[i].yUid, all_ovlp->x[uId].a.a[i].xUid);
if(index == -1) fprintf(stderr, "ERROR 5\n");
all_ovlp->x[all_ovlp->x[uId].a.a[i].yUid].a.a[index].status = FLIP;
}
}
}
}
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
///has been removed as contained
if(purge_g->seq[uId].del || purge_g->seq[uId].c == ALTER_LABLE) continue;
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(all_ovlp->x[uId].a.a[i].status == DELETE) continue;
if(purge_g->seq[all_ovlp->x[uId].a.a[i].xUid].c == ALTER_LABLE||
purge_g->seq[all_ovlp->x[uId].a.a[i].xUid].del||
purge_g->seq[all_ovlp->x[uId].a.a[i].yUid].c == ALTER_LABLE||
purge_g->seq[all_ovlp->x[uId].a.a[i].yUid].del)
{
continue;
}
///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
r = get_hap_arch(&(all_ovlp->x[uId].a.a[i]), ug->u.a[all_ovlp->x[uId].a.a[i].xUid].len,
ug->u.a[all_ovlp->x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
if(r < 0) continue;
p_t = asg_arc_pushp(purge_g);
*p_t = t;
}
}
asg_cleanup(purge_g);
asg_symm(purge_g);
clean_purge_graph(purge_g, keep_rate, 0);
asg_arc_t *av = NULL;
uint32_t n_vtx = purge_g->n_seq*2, nv, w;
for (v = 0; v < n_vtx; ++v)
{
if(purge_g->seq[v>>1].c == ALTER_LABLE || purge_g->seq[v>>1].del) continue;
nv = asg_arc_n(purge_g, v);
av = asg_arc_a(purge_g, v);
if (nv == 0) continue;
for (i = 0; i < nv; ++i)
{
if (av[i].del) continue;
w = av[i].v;
index = get_specific_hap_overlap(&(all_ovlp->x[v>>1]), v>>1, w>>1);
all_ovlp->x[v>>1].a.a[index].status = FLIP;
}
}
p = NULL;
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
for (i = m = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(all_ovlp->x[uId].a.a[i].status != FLIP) continue;
all_ovlp->x[uId].a.a[m] = all_ovlp->x[uId].a.a[i];
if(p == NULL || p->score > all_ovlp->x[uId].a.a[m].score)
{
p = &(all_ovlp->x[uId].a.a[m]);
}
m++;
}
all_ovlp->x[uId].a.n = m;
}
max_score = 0;
if(p && p->score <= 0)
{
max_score = ((p->score)*-1) + 1;
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
all_ovlp->x[uId].a.a[i].score += max_score;
}
}
}
init_contig_phase(all_ovlp, purge_g, drop_ratio);
}
void purge_dups_back(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density,
uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain,
@@ -5243,7 +4989,7 @@ void debug_p_g_t(p_g_t* pg, hap_cov_t *cov, asg_t *read_g)
offset += (uint32_t)u->a[i];
if(i >= sid && i <= eid)
{
if((!!cov->t_ch->is_r_het[u->a[i]>>33]) != t->h_status)
if(cov->t_ch->is_r_het[u->a[i]>>33] != t->h_status)
{
fprintf(stderr, "ERROR-(-3): is_r_het: %u, h_status: %u\n", cov->t_ch->is_r_het[u->a[i]>>33], t->h_status);
}
@@ -5390,15 +5136,15 @@ p_g_t *init_p_g_t_back(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
{
uint32_t v, uId, k, l, offset, l_pos, g_beg_idx, occ, ovlp, tLen, zLen;
uint32_t v, uId, k, l, offset, l_pos, g_beg_idx, occ/**, ovlp, tLen, zLen**/;
p_g_t *pg = NULL; CALLOC(pg, 1);
pg->ug = ug;
asg_t* nsg = pg->ug->g;
ma_utg_t *u = NULL;
p_node_t *t = NULL, *z = NULL;
asg_arc_t *e = NULL;
p_node_t *t = NULL/**, *z = NULL**/;
///asg_arc_t *e = NULL;
p_g_in_t *x = NULL;
pg->pg_het = asg_init();
///pg->pg_het = asg_init();
pg->pg_h_lev = asg_init();
kv_init(pg->pg_het_node);
kv_init(pg->pg_h_lev_idx);
@@ -5417,6 +5163,14 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
pg->pg_h_lev->seq[uId].c = PRIMARY_LABLE;
}
// if(asm_opt.polyploidy <= 2)
// {
// for (v = 0; v < cov->t_ch->r_num; v++)
// {
// if(cov->t_ch->is_r_het[v]&P_HET) cov->t_ch->is_r_het[v] |= S_HET;
// }
// }
for (v = 0; v < nsg->n_seq; v++)
{
uId = v;
@@ -5427,18 +5181,18 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
///fprintf(stderr, "\n+v: %u, pg->pg_het_node.n: %u\n", v, (uint32_t)pg->pg_het_node.n);
for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k)
{
///if(k < u->n && cov->t_ch->is_r_het[u->a[k]>>33] == 0) fprintf(stderr, "sbsbs-uId=%u\n", uId);
if (k == u->n || (!!cov->t_ch->is_r_het[u->a[k]>>33]) != (!!cov->t_ch->is_r_het[u->a[l]>>33]))
///if (k == u->n || (!!cov->t_ch->is_r_het[u->a[k]>>33]) != (!!cov->t_ch->is_r_het[u->a[l]>>33]))
if (k == u->n || cov->t_ch->is_r_het[u->a[k]>>33] != cov->t_ch->is_r_het[u->a[l]>>33])
{
kv_pushp(p_node_t, pg->pg_het_node, &t);
t->c_ug_id = uId;
t->h_status = (!!cov->t_ch->is_r_het[u->a[l]>>33]);
t->h_status = cov->t_ch->is_r_het[u->a[l]>>33];
t->baseBeg = l_pos;
t->baseEnd = offset + read_g->seq[u->a[k-1]>>33].len - 1;
t->nodeBeg = l;
t->nodeEnd = k - 1;
///if(t->b_ug_id == (uint32_t)-1) fprintf(stderr, "xxxx\n");
asg_seq_set(pg->pg_het, pg->pg_het_node.n-1, t->baseEnd+1-t->baseBeg, 0);
///asg_seq_set(pg->pg_het, pg->pg_het_node.n-1, t->baseEnd+1-t->baseBeg, 0);
///fprintf(stderr, "l: %u, k: %u, u->n: %u, t->h_status: %u\n", l, k, u->n, t->h_status);
l = k;
l_pos = offset + (uint32_t)u->a[k-1];
@@ -5451,29 +5205,29 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
x->beg = g_beg_idx; x->occ = occ;
///fprintf(stderr, "-v: %u, pg->pg_het_node.n: %u\n", v, (uint32_t)pg->pg_het_node.n);
if(occ > 1)
{
for (k = g_beg_idx; (k + 1) < pg->pg_het_node.n; ++k)
{
t = &(pg->pg_het_node.a[k]); tLen = t->baseEnd + 1 - t->baseBeg;
z = &(pg->pg_het_node.a[k+1]); zLen = z->baseEnd + 1 - z->baseBeg;
// if(occ > 1)
// {
// for (k = g_beg_idx; (k + 1) < pg->pg_het_node.n; ++k)
// {
// t = &(pg->pg_het_node.a[k]); tLen = t->baseEnd + 1 - t->baseBeg;
// z = &(pg->pg_het_node.a[k+1]); zLen = z->baseEnd + 1 - z->baseBeg;
ovlp = ((MIN(t->baseEnd, z->baseEnd) >= MAX(t->baseBeg, z->baseBeg))?
MIN(t->baseEnd, z->baseEnd) - MAX(t->baseBeg, z->baseBeg) + 1 : 0);
// ovlp = ((MIN(t->baseEnd, z->baseEnd) >= MAX(t->baseBeg, z->baseBeg))?
// MIN(t->baseEnd, z->baseEnd) - MAX(t->baseBeg, z->baseBeg) + 1 : 0);
e = asg_arc_pushp(pg->pg_het);
e->ol = ovlp;
e->ul = (k<<1); e->ul <<= 32; e->ul += (tLen - ovlp);
e->v = ((k+1)<<1); e->del = 0; e->el = e->no_l_indel = e->strong = 1;
// e = asg_arc_pushp(pg->pg_het);
// e->ol = ovlp;
// e->ul = (k<<1); e->ul <<= 32; e->ul += (tLen - ovlp);
// e->v = ((k+1)<<1); e->del = 0; e->el = e->no_l_indel = e->strong = 1;
e = asg_arc_pushp(pg->pg_het);
e->ol = ovlp;
e->ul = ((k+1)<<1)+1; e->ul <<= 32; e->ul += (zLen - ovlp);
e->v = (k<<1)+1; e->del = 0; e->el = e->no_l_indel = e->strong = 1;
}
}
// e = asg_arc_pushp(pg->pg_het);
// e->ol = ovlp;
// e->ul = ((k+1)<<1)+1; e->ul <<= 32; e->ul += (zLen - ovlp);
// e->v = (k<<1)+1; e->del = 0; e->el = e->no_l_indel = e->strong = 1;
// }
// }
}
asg_cleanup(pg->pg_het);
///asg_cleanup(pg->pg_het);
///debug_p_g_t(pg, cov, read_g);
///print_p_g_t_interval(pg, cov);
@@ -5493,6 +5247,147 @@ void destory_p_g_t(p_g_t **pg)
}
}
void partition_contigs(hap_overlaps_list* all_ovlp, ma_ug_t *ug, hap_cov_t *cov, double keep_rate,
int max_hang, int min_ovlp, float drop_ratio, p_g_t *pg)
{
int r, index;
uint32_t v, i, uId, m;
hap_overlaps *p = NULL;
asg_arc_t t, *p_t = NULL;
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
///has been removed as contained
if(pg->pg_h_lev->seq[uId].del || pg->pg_h_lev->seq[uId].c == ALTER_LABLE) continue;
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(all_ovlp->x[uId].a.a[i].type == YCX) continue;
if(all_ovlp->x[uId].a.a[i].type == XCY) continue;
/****************************may have bugs********************************/
if(all_ovlp->x[uId].a.a[i].score <= 0) continue;
/****************************may have bugs********************************/
if(pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].xUid].c == ALTER_LABLE||
pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].xUid].del||
pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].yUid].c == ALTER_LABLE||
pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].yUid].del)
{
continue;
}
///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
r = get_hap_arch(&(all_ovlp->x[uId].a.a[i]), ug->u.a[all_ovlp->x[uId].a.a[i].xUid].len,
ug->u.a[all_ovlp->x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
if(r < 0) continue;
p_t = asg_arc_pushp(pg->pg_h_lev);
*p_t = t;
}
}
asg_cleanup(pg->pg_h_lev);
asg_symm(pg->pg_h_lev);
clean_purge_graph(pg->pg_h_lev, keep_rate, 0);
asg_arc_t *av = NULL;
uint32_t n_vtx = pg->pg_h_lev->n_seq<<1, nv, a, b;
for (v = 0; v < n_vtx; v++)
{
av = asg_arc_a(pg->pg_h_lev, v);
nv = asg_arc_n(pg->pg_h_lev, v);
for (i = 0; i < nv; ++i)
{
if(av[i].del) continue;
a = av[i].ul>>33;
b = av[i].v>>1;
index = get_specific_hap_overlap(&(all_ovlp->x[a]), a, b);
p = &(all_ovlp->x[a].a.a[index]);
p->status = MIXED;
index = get_specific_hap_overlap(&(all_ovlp->x[b]), b, a);
p = &(all_ovlp->x[b].a.a[index]);
p->status = MIXED;
}
}
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
for (i = m = 0; i < all_ovlp->x[uId].a.n; i++)
{
p = (&all_ovlp->x[uId].a.a[i]);
/****************************may have bugs********************************/
if(p->score <= 0) continue;
/****************************may have bugs********************************/
if(p->type == YCX || p->type == XCY || p->status == MIXED)
{
all_ovlp->x[uId].a.a[m] = (*p);
all_ovlp->x[uId].a.a[m].status = SELF_EXIST;
m++;
}
}
all_ovlp->x[uId].a.n = m;
}
for (v = 0; v < all_ovlp->num; v++)
{
uId = v; p = NULL;
if(all_ovlp->x[uId].a.n == 0) continue;
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(p == NULL || p->score < all_ovlp->x[uId].a.a[i].score)
{
p = &(all_ovlp->x[uId].a.a[i]);
}
}
if(!p) continue;
for (i = m = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(all_ovlp->x[uId].a.a[i].type == YCX || all_ovlp->x[uId].a.a[i].type == XCY)
{
if(filter_secondary_chain(p->score, all_ovlp->x[uId].a.a[i].score, keep_rate))
{
all_ovlp->x[uId].a.a[m] = all_ovlp->x[uId].a.a[i];
m++;
}
}
else
{
all_ovlp->x[uId].a.a[m] = all_ovlp->x[uId].a.a[i];
m++;
}
}
all_ovlp->x[uId].a.n = m;
}
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
for (i = m = 0; i < all_ovlp->x[uId].a.n; i++)
{
p = (&all_ovlp->x[uId].a.a[i]);
index = get_specific_hap_overlap(&(all_ovlp->x[p->yUid]), p->yUid, p->xUid);
if(index != -1)
{
all_ovlp->x[uId].a.a[m] = (*p);
m++;
}
}
all_ovlp->x[uId].a.n = m;
}
}
void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density,
uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain,
@@ -5573,6 +5468,13 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);
filter_hap_overlaps_by_length(&all_ovlp, purege_minLen);
if(collect_p_trans)
{
pt_solve(&all_ovlp, cov->t_ch, ug, read_g, 0.8, R_INF.trio_flag);
goto end_coverage;
}
pg = init_p_g_t(ug, cov, read_g);
///normalize_hap_overlaps(&all_ovlp, &back_all_ovlp);
Purge_Dups.h
+28 -1
View File
@@ -13,6 +13,33 @@
#define REAL_ALTER_THRES 0.1
#define CHAIN_FILTER_RATE 0.7
typedef struct {
uint8_t rev;
uint8_t type;
uint8_t status;
uint32_t x_beg_pos;
uint32_t x_end_pos;
uint32_t y_beg_pos;
uint32_t y_end_pos;
uint32_t x_beg_id;
uint32_t x_end_id;
uint32_t y_beg_id;
uint32_t y_end_id;
uint32_t xUid;
uint32_t yUid;
uint32_t weight;
long long score;
}hap_overlaps;
typedef struct {
kvec_t(hap_overlaps) a;
}kvec_hap_overlaps;
typedef struct {
kvec_hap_overlaps* x;
uint32_t num;
}hap_overlaps_list;
void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density,
uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain,
@@ -25,6 +52,6 @@ hap_cov_t* init_hap_cov_t(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* sources, R
ma_sub_t *coverage_cut, int max_hang, int min_ovlp, uint32_t is_collect_trans);
void destory_hap_cov_t(hap_cov_t **x);
void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd);
int get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn);
#endif
hic.cpp
+225 -13
View File
@@ -19,7 +19,7 @@ KDQ_INIT(uint64_t)
#define OFFSET_RATE 0.000000001
#define OFFSET_SECOND_RATE 0.0000000001
#define SCALL 10000
#define OFFSET_RATE_MAX_W 20.8286263517*SCALL
#define OFFSET_RATE_MAX_W 6.90675477865*SCALL
#define OFFSET_RATE_MIN_W 4.0000003e-10*SCALL
#define HIC_COUNTER_BITS 12
@@ -4143,8 +4143,8 @@ void print_hc_links(hc_links* link, int dir, H_partition* hap)
}
}
fprintf(stderr, "self-utg%.6dl\tFW:%f\tRW:%f\tRT:%f\n**************************************************\n",
(int)(i+1), f_w, r_w, (f_w+r_w) != 0? r_w/(f_w+r_w):0);
// fprintf(stderr, "self-utg%.6dl\tFW:%f\tRW:%f\tRT:%f\n**************************************************\n",
// (int)(i+1), f_w, r_w, (f_w+r_w) != 0? r_w/(f_w+r_w):0);
}
}
@@ -7437,6 +7437,7 @@ void update_bsg(asg_t *bsg, kvec_asg_arc_t_warp* edges)
}
void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub)
{
double index_time = yak_realtime();
ma_ug_t *bub_ug = bub->b_ug;
asg_t *bsg = bub->b_g;
uint32_t k, i, v, n_vx, new_bub;
@@ -7605,11 +7606,13 @@ void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub)
kv_destroy(res_btg.a); kv_destroy(res_utg.a); kv_destroy(edges.a);
///print_debug_bubble_graph(bub, ug, asm_opt.output_file_name);
fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
}
void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint32_t is_end)
{
double index_time = yak_realtime();
if(bub->b_ug) ma_ug_destroy(bub->b_ug);
if(bub->chain_weight.a) kv_destroy(bub->chain_weight);
kvec_t_u32_warp broken;
@@ -7828,7 +7831,7 @@ void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint
}
}
**/
fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
}
void set_b_utg_weight_flag(bubble_type* bub, buf_t* b, uint32_t v, uint8_t* vis_flag, uint32_t flag, uint32_t* occ)
@@ -8187,6 +8190,7 @@ int double_check_bub_branch(asg_arc_t *t, ma_ug_t *bs_ug, double *e_w, uint32_t
void clean_bubble_chain_by_HiC(ma_ug_t* ug, hc_links* link, bubble_type* bub)
{
double index_time = yak_realtime();
ma_ug_t *bs_ug = bub->b_ug;
uint32_t v, u, i, m, max_i, nv, rv, n_vx, root, flag_pri = 1, flag_aux = 2, flag_ava = 4, occ;
double w, cutoff = 2/**, max_w_cutoff = MAX(MIN(100*OFFSET_RATE_MIN_W, OFFSET_RATE_MAX_W/100), OFFSET_RATE_MIN_W)**/;
@@ -8383,11 +8387,13 @@ void clean_bubble_chain_by_HiC(ma_ug_t* ug, hc_links* link, bubble_type* bub)
free(vis); free(is_vis); free(is_used); free(dedup); free(b.b.a); free(e_w); free(e_occ);
kv_destroy(stack.a); kv_destroy(result.a); kv_destroy(res_utg.a); kv_destroy(edges.a);
ma_ug_destroy(back_bs_ug);
fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
}
void append_boundary_chain(ma_ug_t* ug, hc_links* link, bubble_type* bub)
{
double index_time = yak_realtime();
ma_ug_t *bs_ug = bub->b_ug;
uint32_t v, u, i, k, beg_idx, m, nv, n_vx, flag_pri = 1, flag_aux = 2, flag_ava = 4;
uint32_t root, rId_0, ori_0, root_0, new_bub;
@@ -8538,6 +8544,7 @@ void append_boundary_chain(ma_ug_t* ug, hc_links* link, bubble_type* bub)
// }
// }
/*******************************for debug************************************/
fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
}
int cmp_chain_hic_w_weight(const void * a, const void * b)
@@ -9408,12 +9415,12 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link,
return 0;
}
// for (i = 0; i < dis->n; i++)
// {
// if(i > 0 && dis->a[i].beg != dis->a[i-1].end) fprintf(stderr, "ERROR: dis->a[i].beg: %lu, dis->a[i-1].end: %lu\n", dis->a[i].beg, dis->a[i-1].end);
// fprintf(stderr, "beg: %lu, end: %lu, cnt_0: %lu, cnt_1: %lu, error_rate: %f\n",
// dis->a[i].beg, dis->a[i].end, dis->a[i].cnt_0, dis->a[i].cnt_1, (double)(dis->a[i].cnt_1)/(double)(dis->a[i].cnt_1 + dis->a[i].cnt_0));
// }
for (i = 0; i < dis->n; i++)
{
if(i > 0 && dis->a[i].beg != dis->a[i-1].end) fprintf(stderr, "ERROR: dis->a[i].beg: %lu, dis->a[i-1].end: %lu\n", dis->a[i].beg, dis->a[i-1].end);
fprintf(stderr, "beg: %lu, end: %lu, cnt_0: %lu, cnt_1: %lu, error_rate: %f\n",
dis->a[i].beg, dis->a[i].end, dis->a[i].cnt_0, dis->a[i].cnt_1, (double)(dis->a[i].cnt_1)/(double)(dis->a[i].cnt_1 + dis->a[i].cnt_0));
}
LeastSquare_advance(dis, idx, med);
// fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n",
@@ -9441,6 +9448,196 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link,
}
int get_trans_rate_function_advance(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, MT* M, H_partition* hap, trans_idx* dis)
{
kvec_t(uint64_t) buf;
kv_init(buf);
uint64_t beg, end, cnt[2];
uint64_t k, i, t_d, r_idx, f_idx, med = (uint64_t)-1;
int beg_status, end_status;
for (i = 0; i < link->a.n; i++)
{
for (k = 0; k < link->a.a[i].e.n; k++)
{
link->a.a[i].e.a[k].dis = (uint64_t)-1;
}
}
fill_utg_distance_multi(idx, link, M, bub);
buf.n = 0;
for (k = 0; k < hits->a.n; ++k)
{
beg = ((hits->a.a[k].s<<1)>>(64 - idx->uID_bits));
end = ((hits->a.a[k].e<<1)>>(64 - idx->uID_bits));
if(IF_HOM(beg, *bub)) continue;
if(IF_HOM(end, *bub)) continue;
t_d = get_hic_distance(&(hits->a.a[k]), link, idx);
if(t_d == (uint64_t)-1) continue;
if(beg == end)
{
t_d = (t_d << 1);
}
else
{
beg_status = get_phase_status(hap, beg);
if(beg_status != 1 && beg_status != -1) continue;
end_status = get_phase_status(hap, end);
if(end_status != 1 && end_status != -1) continue;
if(beg_status != end_status)
{
t_d = (t_d << 1) + 1;
}
else
{
t_d = (t_d << 1);
}
}
kv_push(uint64_t, buf, t_d);
}
///might have bias, we may not use right linkage larger than trans rc linkage
radix_sort_hc64(buf.a, buf.a+buf.n);
for (k = 0, r_idx = f_idx = (uint64_t)-1; k < buf.n; k++)
{
if((buf.a[k]&1) == 0) r_idx = k;
if((buf.a[k]&1) == 1) f_idx = k;
}
buf.n = MIN(r_idx, f_idx);
trans_p_t* p = NULL;
dis->n = 0;
uint64_t bin_size = MIN(2250, buf.n>>8), m;
i = 0;
while (i < buf.n)
{
kv_pushp(trans_p_t, *dis, &p);
p->beg = i;
cnt[0] = cnt[1] = 0;
k = MIN(i+bin_size, buf.n);
for (; i < k; i++)
{
cnt[buf.a[i]&1]++;
}
p->end = i;
p->cnt_0 = cnt[0];
p->cnt_1 = cnt[1];
}
i = m = 0; med = (uint64_t)-1;
while(i < dis->n)
{
if(dis->a[i].cnt_0 > 0 && dis->a[i].cnt_1 > 0)
{
dis->a[m] = dis->a[i];
m++;
i++;
continue;
}
if(med == (uint64_t)-1) med = buf.a[dis->a[i].beg]>>1;
k = i; cnt[0] = cnt[1] = 0;
for (; i < dis->n; i++)
{
cnt[0] += dis->a[i].cnt_0;
cnt[1] += dis->a[i].cnt_1;
if(cnt[0] > 0 && cnt[1] > 0) break;
}
if(i < dis->n)
{
dis->a[m].cnt_0 = cnt[0];
dis->a[m].cnt_1 = cnt[1];
dis->a[m].beg = dis->a[k].beg;
dis->a[m].end = dis->a[i].end;
m++;
i++;
continue;
}
cnt[0] -= dis->a[k].cnt_0;
cnt[1] -= dis->a[k].cnt_1;
while (1)
{
cnt[0] += dis->a[k].cnt_0;
cnt[1] += dis->a[k].cnt_1;
if(cnt[0] > 0 && cnt[1] > 0) break;
if(k == 0)
{
k = (uint64_t)-1;
break;
}
k--;
}
if(k != (uint64_t)-1)
{
dis->a[m].cnt_0 = cnt[0];
dis->a[m].cnt_1 = cnt[1];
dis->a[m].end = dis->a[i-1].end;
continue;
}
m = 0;
break;
}
dis->n = m;
if(dis->n == 0)
{
kv_destroy(buf);
return 0;
}
for (i = 0; i < dis->n; i++)
{
dis->a[i].beg = buf.a[dis->a[i].beg]>>1;
dis->a[i].end = (buf.a[dis->a[i].end-1]>>1);
}
for (i = 0; i < dis->n - 1; i++)
{
dis->a[i].end += ((dis->a[i+1].beg - dis->a[i].end)/2);
dis->a[i+1].beg = dis->a[i].end;
}
for (i = 0; i < dis->n; i++)
{
if(i > 0 && dis->a[i].beg != dis->a[i-1].end) fprintf(stderr, "ERROR: dis->a[i].beg: %lu, dis->a[i-1].end: %lu\n", dis->a[i].beg, dis->a[i-1].end);
fprintf(stderr, "beg: %lu, end: %lu, cnt_0: %lu, cnt_1: %lu, error_rate: %f\n",
dis->a[i].beg, dis->a[i].end, dis->a[i].cnt_0, dis->a[i].cnt_1, (double)(dis->a[i].cnt_1)/(double)(dis->a[i].cnt_1 + dis->a[i].cnt_0));
}
LeastSquare_advance(dis, idx, med);
// fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n",
// (double)idx->a, (double)idx->b, (double)idx->frac, med);
dis->max = dis->a[dis->n-1].end;
kv_destroy(buf);
if(idx->a < 0) idx->a = 0;
if(idx->a == 0)
{
idx->b = MAX((((double)(dis->a[dis->n-1].cnt_1))/((double)(dis->a[dis->n-1].cnt_0 + dis->a[dis->n-1].cnt_1))), idx->b);
}
if(idx->b < 0 && get_trans(idx, dis->max) < 0)
{
idx->b = ((double)(dis->a[dis->n-1].cnt_1))/((double)(dis->a[dis->n-1].cnt_0 + dis->a[dis->n-1].cnt_1));
}
// fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n",
// (double)idx->a, (double)idx->b, (double)idx->frac, med);
return 1;
}
void init_hic_p(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub,
kvec_hc_edge* back_hc_edge, MT* M, H_partition* hap, uint32_t ignore_dis)
{
@@ -9450,7 +9647,7 @@ kvec_hc_edge* back_hc_edge, MT* M, H_partition* hap, uint32_t ignore_dis)
if(bub->round_id > 0 && ignore_dis == 0)
{
is_comples_weight = get_trans_rate_function(idx, hits, link, bub, M, hap, &dis);
is_comples_weight = get_trans_rate_function_advance(idx, hits, link, bub, M, hap, &dis);
}
@@ -10304,6 +10501,7 @@ void print_contig_partition(H_partition* hap, const char* debug)
void adjust_contig_partition(H_partition* hap, hc_links* link)
{
double index_time = yak_realtime();
uint32_t i, k, *h0, h0_n, *h1, h1_n;
uint32_t h0_status[4], h1_status[4], h0_status_max;
int h0_h, h1_h;
@@ -10404,6 +10602,7 @@ void adjust_contig_partition(H_partition* hap, hc_links* link)
hap->g_p.a[i].weight[1] = get_cluster_weight(hap, link, h1, h1_n);
hap->g_p.a[i].weight_convex = get_cluster_inner_weight(hap, link, h0, h0_n, h1, h1_n);
}
fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
}
uint32_t get_weightest_uid(uint32_t* a, uint32_t n, H_partition* hap, uint8_t* hap_label_flag,
@@ -11733,6 +11932,7 @@ int graph_bipartiteness(uint32_t* b_a, uint32_t b_a_n, uint8_t *color, hc_links*
void assign_per_unitig_G_partition(G_partition* g_p, uint64_t hap_n, hc_links* link, bubble_type* bub,
uint32_t bubble_first)
{
double index_time = yak_realtime();
reset_G_partition(g_p, hap_n);
partition_warp* res = NULL;
@@ -11872,7 +12072,7 @@ uint32_t bubble_first)
g_p->index[res->a.a[k]] = (g_p->index[res->a.a[k]] << 1) + 1;
}
}
fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
}
typedef struct {
@@ -12093,6 +12293,7 @@ bub_sort_vec* w_stack, uint8_t* vis, uint32_t* set_hap, uint32_t n_utg, uint32_t
void sort_bubble_ele(G_partition* g_p, hc_links* link, bubble_type* bub, uint32_t n_utg)
{
double index_time = yak_realtime();
kvec_t_u64_warp stack; kv_init(stack.a);
bub_sort_vec w_stack; kv_init(w_stack);
uint8_t* vis = NULL; MALLOC(vis, n_utg);
@@ -12106,6 +12307,7 @@ void sort_bubble_ele(G_partition* g_p, hc_links* link, bubble_type* bub, uint32_
}
kv_destroy(stack.a); kv_destroy(w_stack); free(vis); free(set_hap);
fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
}
uint32_t init_contig_partition(H_partition* hap, ha_ug_index* idx, bubble_type* bub, hc_links* link)
@@ -13214,22 +13416,32 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx)
bub.round_id = 0; bub.n_round = 2;
for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++)
{
fprintf(stderr, "0********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime());
identify_bubbles(idx->ug, &bub, idx->cov->t_ch->is_r_het);
fprintf(stderr, "1********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime());
if(bub.round_id == 0)
{
fprintf(stderr, "2********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime());
collect_hc_links(sl.idx, &sl.hits, &link, &bub, &M);
fprintf(stderr, "3********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime());
collect_hc_reverse_links(&link, idx->ug, &bub);
fprintf(stderr, "4********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime());
}
init_hic_p((ha_ug_index*)sl.idx, &sl.hits, &link, &bub, &back_hc_edge, &M, &hap, 0);
fprintf(stderr, "5********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime());
///init_hic_p_new((ha_ug_index*)sl.idx, &sl.hits, idx->link, &bub, &back_hc_edge, &M);
reset_H_partition(&hap, (bub.round_id == 0? 1 : 0));
fprintf(stderr, "6********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime());
init_contig_partition(&hap, idx, &bub, &link);
fprintf(stderr, "7********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime());
phasing_improvement(&hap, &(hap.g_p), idx, &bub, &link);
fprintf(stderr, "8********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime());
label_unitigs(&(hap.g_p), idx->ug);
fprintf(stderr, "9********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime());
///print_hc_links(idx->link, 0, &hap);
}
///print_hc_links(idx->link, 0, &hap);
print_hc_links(&link, 0, &hap);
cluster_contigs(&bub, idx, &sl.hits, &M, &hap, &link);
partig.cpp
+732
View File
@@ -0,0 +1,732 @@
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include <stdlib.h>
#include "partig.h"
#include "Purge_Dups.h"
#include "Correct.h"
#include "ksort.h"
#define generic_key(x) (x)
KRADIX_SORT_INIT(gfa64, uint64_t, generic_key, 8)
#define pt_a(x, id) ((x).ma.a + ((x).idx.a[(id)]>>32))
#define pt_n(x, id) ((uint32_t)((x).idx.a[(id)]))
typedef struct {
int32_t topn;
int32_t n_perturb;
uint64_t seed;
double f_perturb;
} pt_svopt_t;
typedef struct {
///uint32_t m, n, *shuffled;
uint32_t *shuffled;
uint32_t off, size; // offset in pt_match_t::cc; size of the component
uint64_t *buf;
kvec_t(uint64_t) a;
int8_t *s, *s_tmp;///s is the status (haplotype) of each unitig: for backup
} solve_aux_t;
static inline uint64_t kr_splitmix64(uint64_t x)
{
uint64_t z = (x += 0x9E3779B97F4A7C15ULL);
z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL;
z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL;
return z ^ (z >> 31);
}
static inline double kr_drand_r(uint64_t *x)
{
union { uint64_t i; double d; } u;
*x = kr_splitmix64(*x);
u.i = 0x3FFULL << 52 | (*x) >> 12;
return u.d - 1.0;
}
void pt_svopt_init(pt_svopt_t *opt)
{
memset(opt, 0, sizeof(pt_svopt_t));
opt->seed = 11;
opt->topn = 1<<30;
opt->n_perturb = 1000;
opt->f_perturb = 0.1;
}
static void pt_pdist_idx(pt_match_t *ma)
{
uint32_t st, i;
kv_resize(uint64_t, ma->idx, ma->n_seq);
ma->idx.n = ma->n_seq;
memset(ma->idx.a, 0, ma->idx.n*sizeof(uint64_t));
for (st = 0, i = 1; i <= ma->ma.n; ++i)
if (i == ma->ma.n || ma->ma.a[i].sid[0] != ma->ma.a[st].sid[0])
ma->idx.a[ma->ma.a[st].sid[0]] = (uint64_t)st << 32 | (i - st), st = i;
}
static pt_match1_t *pt_pdist(const pt_match_t *ma, uint32_t sid1, uint32_t sid2)
{
pt_match1_t *o = pt_a(*ma, sid1);
uint32_t n = pt_n(*ma, sid1), k;
for (k = 0; k < n; ++k)
if (o[k].sid[1] == sid2)
return &(o[k]);
return NULL;
}
static void normalize_pdist(pt_match1_t *a, pt_match1_t *b)
{
if(a->w >= b->w)
{
b->sid[0] = a->sid[1];
b->sid[1] = a->sid[0];
b->w = a->w;
}
else
{
a->sid[0] = b->sid[1];
a->sid[1] = b->sid[0];
a->w = b->w;
}
}
uint32_t pt_pdist_symm(pt_match_t *ma)
{
uint8_t *del = NULL;
uint32_t i, k, n = 0;
pt_match1_t *t = NULL, *m = NULL;
CALLOC(del, ma->ma.n);
for (i = 0; i < ma->ma.n; ++i) {
m = &ma->ma.a[i];
if (m->sid[0] == m->sid[1])
{
del[i] = 1, ++n;///self overlap
continue;
}
t = pt_pdist(ma, m->sid[1], m->sid[0]);
if(!t)
{
del[i] = 1, ++n;///self overlap
continue;
}
normalize_pdist(m, t);
}
if (n > 0) {
for (i = k = 0; i < ma->ma.n; ++i)
if (!del[i]) ma->ma.a[k++] = ma->ma.a[i];
ma->ma.n = k;
pt_pdist_idx(ma);
}
free(del);
return n;
}
static void pt_pdist_flt(pt_match_t *ma, uint32_t min_cnt, double drop_thres)
{
uint32_t i, k, j, n, max, o;
uint8_t *mark = NULL;
CALLOC(mark, ma->ma.n);
for (i = 0; i < ma->n_seq; ++i) {
o = ma->idx.a[i] >> 32;
n = (uint32_t)ma->idx.a[i];
max = 0;
if (n == 0) continue;
for (j = o; j < o + n; ++j)
max = max > ma->ma.a[j].w? max : ma->ma.a[j].w;
for (j = o; j < o + n; ++j)
if (ma->ma.a[j].w >= (max*drop_thres) || ma->ma.a[j].w + min_cnt >= max)
mark[j] = 1;
}
for (i = 0; i < ma->ma.n; ++i)
{
if (mark[i] == 0) continue;
o = ma->idx.a[ma->ma.a[i].sid[1]]>>32;
n = (uint32_t)ma->idx.a[ma->ma.a[i].sid[1]];
for (j = o; j < o + n; ++j)
{
if (ma->ma.a[j].sid[1] == ma->ma.a[i].sid[0]) mark[j] = 1;
}
}
for (i = k = 0; i < ma->ma.n; ++i)
if (mark[i]) ma->ma.a[k++] = ma->ma.a[i];
ma->ma.n = k;
free(mark);
pt_pdist_idx(ma);
pt_pdist_symm(ma);
}
pt_match_t *init_pt_match_t(hap_overlaps_list* ha, pt_g_t *x, double f_rate)
{
pt_match_t *p = NULL; CALLOC(p, 1); p->n_seq = x->ug->g->n_seq;
kv_init(p->idx); kv_init(p->ma);
uint32_t v, i, k, qn, tn, qs, qe, ts, te, occ, as, ae;
uint64_t hetLen, homLen, oLen;
pt_node_t *a = NULL;
pt_match1_t *ma = NULL;
for (v = 0; v < ha->num; v++)
{
for (i = 0; i < ha->x[v].a.n; i++)
{
if(ha->x[v].a.a[i].score <= 0) continue;
if(ha->x[v].a.a[i].xUid == ha->x[v].a.a[i].yUid) continue;
/*****************qn*****************/
qn = ha->x[v].a.a[i].xUid;
qs = ha->x[v].a.a[i].x_beg_pos;
qe = ha->x[v].a.a[i].x_end_pos - 1;
a = x->p.a + x->p_idx.a[qn];
occ = x->p_idx.a[qn+1] - x->p_idx.a[qn];
for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
{
as = a[k].bS;
ae = a[k].bE;
oLen = ((MIN(qe, ae) >= MAX(qs, as))? MIN(qe, ae) - MAX(qs, as) + 1 : 0);
if(homLen + hetLen > 0 && oLen == 0) break;
if(oLen == 0) continue;
if(a[k].hs == N_HET)
{
homLen += oLen;
}
else if(asm_opt.polyploidy <= 2 && (a[k].hs&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].hs&S_HET))
{
homLen += oLen;
}
else
{
hetLen += oLen;
}
}
if(hetLen <= ((hetLen + homLen)*f_rate)) continue;
/*****************qn*****************/
/*****************tn*****************/
tn = ha->x[v].a.a[i].yUid;
ts = ha->x[v].a.a[i].y_beg_pos;
te = ha->x[v].a.a[i].y_end_pos - 1;
a = x->p.a + x->p_idx.a[tn];
occ = x->p_idx.a[tn+1] - x->p_idx.a[tn];
for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
{
as = a[k].bS;
ae = a[k].bE;
oLen = ((MIN(te, ae) >= MAX(ts, as))? MIN(te, ae) - MAX(ts, as) + 1 : 0);
if(homLen + hetLen > 0 && oLen == 0) break;
if(oLen == 0) continue;
if(a[k].hs == N_HET)
{
homLen += oLen;
}
else if(asm_opt.polyploidy <= 2 && (a[k].hs&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].hs&S_HET))
{
homLen += oLen;
}
else
{
hetLen += oLen;
}
}
if(hetLen <= ((hetLen + homLen)*f_rate)) continue;
/*****************tn*****************/
kv_pushp(pt_match1_t, p->ma, &ma);
ma->sid[0] = ha->x[v].a.a[i].xUid;
ma->sid[1] = ha->x[v].a.a[i].yUid;
ma->w = ha->x[v].a.a[i].score;
}
}
pt_pdist_idx(p);
pt_pdist_symm(p);
pt_pdist_flt(p, 5, 0.5);
return p;
}
void debug_pt_g_t(pt_g_t *pg)
{
fprintf(stderr, "0----------[M::%s]----------\n", __func__);
uint32_t i, offset, v, sid, eid, spos, epos, p_status, p_uid, occ;
ma_utg_t *u = NULL;
pt_node_t *a = NULL, *t = NULL;
for (v = 0; v < pg->ug->u.n; v++)
{
a = pg->p.a + pg->p_idx.a[v];
occ = pg->p_idx.a[v+1] - pg->p_idx.a[v];
for (i = 0; i < occ; i++)
{
if(a[i].uID != v) fprintf(stderr, "ERROR-s\n");
}
}
for (v = 0, p_status = (uint32_t)-1, p_uid = (uint32_t)-1; v < pg->p.n; v++)
{
t = &(pg->p.a[v]);
sid = t->nS;
eid = t->nE;
spos = t->bS;
epos = t->bE;
if(p_uid == t->uID && p_status == t->hs)
{
fprintf(stderr, "ERROR-a\n");
}
p_status = t->hs;
p_uid = t->uID;
u = &(pg->ug->u.a[t->uID]);
for (i = offset = 0; i < u->n; i++)
{
if(i == sid)
{
if(spos != offset)
{
fprintf(stderr, "ERROR-b\n");
}
}
if(i == eid)
{
if(epos != (offset+pg->rg->seq[u->a[i]>>33].len - 1))
{
fprintf(stderr, "ERROR-c, real end: %u\n",
(uint32_t)(offset+pg->rg->seq[u->a[i]>>33].len - 1));
}
}
offset += (uint32_t)u->a[i];
if(i >= sid && i <= eid)
{
if(pg->t_ch->is_r_het[u->a[i]>>33] != t->hs)
{
fprintf(stderr, "ERROR-d: is_r_het: %u, h_status: %u\n", pg->t_ch->is_r_het[u->a[i]>>33], t->hs);
}
}
}
}
pt_match1_t *o = NULL, *s = NULL;
uint32_t k, n, cnt;
for (i = 0; i < pg->e->n_seq; ++i)
{
o = pt_a(*(pg->e), i); n = pt_n(*(pg->e), i);
for (k = 0; k < n; ++k)
{
if(o[k].sid[0] != i) fprintf(stderr, "ERROR-g\n");
s = pt_pdist(pg->e, o[k].sid[1], o[k].sid[0]);
if(!s) fprintf(stderr, "ERROR-e\n");
if(s)
{
if(!(s->sid[0] == o[k].sid[1] && s->sid[1] == o[k].sid[0] && s->w == o[k].w))
{
fprintf(stderr, "ERROR-f\n");
}
}
}
for (k = cnt = 0; k < pg->e->ma.n; ++k)
{
if(pg->e->ma.a[k].sid[0] == i) cnt++;
}
if(cnt != n) fprintf(stderr, "ERROR-h\n");
}
}
pt_g_t *init_pt_g_t(hap_overlaps_list* ovlp, trans_chain* t_ch, ma_ug_t *ug, asg_t *read_g, double f_rate)
{
uint32_t v, l, k, offset, l_pos;
pt_g_t *p = NULL; CALLOC(p, 1);
asg_t* nsg = ug->g;
ma_utg_t *u = NULL;
pt_node_t *t = NULL;
p->ug = ug;
p->rg = read_g;
p->t_ch = t_ch;
kv_init(p->info); p->info.n = p->info.m = p->ug->g->n_seq; CALLOC(p->info.a, p->info.n);
kv_init(p->p);
kv_init(p->p_idx); kv_push(uint32_t, p->p_idx, 0);
for (v = 0; v < nsg->n_seq; v++)
{
if(nsg->seq[v].del || nsg->seq[v].c == ALTER_LABLE) continue;
u = &(ug->u.a[v]);
for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k)
{
if (k == u->n || t_ch->is_r_het[u->a[k]>>33] != t_ch->is_r_het[u->a[l]>>33])
{
kv_pushp(pt_node_t, p->p, &t);
t->uID = v;
t->hs = t_ch->is_r_het[u->a[l]>>33];
t->bS = l_pos;
t->bE = offset + read_g->seq[u->a[k-1]>>33].len - 1;
t->nS = l;
t->nE = k - 1;
l = k;
l_pos = offset + (uint32_t)u->a[k-1];
}
offset += (uint32_t)u->a[k-1];
}
kv_push(uint32_t, p->p_idx, p->p.n);
}
p->e = init_pt_match_t(ovlp, p, f_rate);
return p;
}
void destory_pt_g_t(pt_g_t **p)
{
if(!p || !(*p)) return;
kv_destroy((*p)->p);
kv_destroy((*p)->info);
kv_destroy((*p)->p_idx);
kv_destroy((*p)->e->idx);
kv_destroy((*p)->e->ma);
free((*p)->e->cc);
free((*p)->e);
free((*p));
}
uint64_t *pt_cc_core(const pt_match_t *ma)
{
uint32_t i, x, y, *flag;
uint64_t *group;
pt_match1_t *o = NULL;
kvec_t(uint32_t) stack; kv_init(stack);
MALLOC(flag, ma->n_seq);
for (i = 0; i < ma->n_seq; ++i)
flag[i] = (uint32_t)-1;
// connected componets
for (i = 0; i < ma->n_seq; ++i) {
if (flag[i] != (uint32_t)-1) continue;
stack.n = 0;
kv_push(uint32_t, stack, i);
while (stack.n > 0) {
uint32_t k, j, n;
stack.n--;
k = stack.a[stack.n];
flag[k] = i;///group id
// n = (uint32_t)ma->idx[k];
// s = ma->idx[k] >> 32;
o = pt_a(*ma, k);
n = pt_n(*ma, k);
for (j = 0; j < n; ++j) {
uint32_t t = o[j].sid[1];
if (flag[t] != (uint32_t)-1) continue;
// if (ns == ms) PT_EXPAND(stack, ms);
// stack[ns++] = t;
kv_push(uint32_t, stack, t);
}
}
}
kv_destroy(stack);
// precalculate the size of each group
CALLOC(group, ma->n_seq);
for (i = 0; i < ma->n_seq; ++i)
group[i] = (uint64_t)flag[i] << 32 | i;
radix_sort_gfa64(group, group + ma->n_seq);
for (i = 1, x = y = 0; i <= ma->n_seq; ++i) {
if (i == ma->n_seq || group[i]>>32 != group[x]>>32) {
uint32_t j;
for (j = x; j < i; ++j)
group[j] = (uint64_t)y << 32 | (uint32_t)group[j];///(group id)|first element in this group
++y, x = i;
}
}
free(flag);
return group;
}
void pt_cc(pt_match_t *ma)
{
ma->cc = pt_cc_core(ma);
}
///inspect top INT edges
static int64_t pt_score(const pt_match_t *ma, uint32_t topn, solve_aux_t *aux)
{
uint32_t i;
int64_t z = 0;
for (i = 0; i < aux->size; ++i) {///aux->size: how many unitigs in this group
uint32_t k = (uint32_t)ma->cc[aux->off + i];
uint32_t o = ma->idx.a[k] >> 32;
uint32_t n = (uint32_t)ma->idx.a[k], j;
for (j = 0; j < n; ++j)
aux->buf[j] = (uint64_t)((uint32_t)-1 - ma->ma.a[o + j].w) << 32 | (o + j);
radix_sort_gfa64(aux->buf, aux->buf + n);
for (j = 0; j < n && j < topn; ++j) {
const pt_match1_t *m = &ma->ma.a[(uint32_t)aux->buf[j]];
z += -(int64_t)m->w * aux->s[m->sid[0]] * aux->s[m->sid[1]];
}
}
return z;
}
static int64_t pt_solve1_init_phase(const pt_match_t *ma, int32_t topn, uint64_t *x, solve_aux_t *aux)
{
uint32_t i;
aux->a.n = 0;
for (i = 0; i < aux->size; ++i) {
uint32_t k = (uint32_t)ma->cc[aux->off + i];///unitig id
uint32_t o = ma->idx.a[k] >> 32;///group id
uint32_t n = (uint32_t)ma->idx.a[k], j;
aux->shuffled[i] = k;///init shuffled status
for (j = 0; j < n; ++j) {
const pt_match1_t *m = &ma->ma.a[o + j];
///first is weight, second part is edge id
kv_push(uint64_t, aux->a, (uint64_t)((uint32_t)-1 - m->w) << 32 | (o + j));
}
}
radix_sort_gfa64(aux->a.a, aux->a.a + aux->a.n);///sort all edges in this group
///randomly assign each unitig
for (i = 0; i < aux->a.n; ++i) { // from the strongest edge to the weakest
const pt_match1_t *m = &ma->ma.a[(uint32_t)aux->a.a[i]];
///sid[0]: query id, sid[1]: target id
///the initial results of aux->s is 0
if (aux->s[m->sid[0]] == 0 && aux->s[m->sid[1]] == 0) {
*x = kr_splitmix64(*x);// random number
aux->s[m->sid[0]] = *x&1? 1 : -1;
aux->s[m->sid[1]] = -aux->s[m->sid[0]];
} else if (aux->s[m->sid[0]] == 0) {
aux->s[m->sid[0]] = -aux->s[m->sid[1]];
} else if (aux->s[m->sid[1]] == 0) {
aux->s[m->sid[1]] = -aux->s[m->sid[0]];
}
}
return pt_score(ma, topn, aux);
}
///size is how many unitigs in this group
static void ks_shuffle_uint32_t(size_t n, uint32_t a[], uint64_t *x)
{
size_t i, j;
for (i = n; i > 1; --i) {
uint32_t tmp;
j = (size_t)(kr_drand_r(x) * i);///semms
tmp = a[j]; a[j] = a[i-1]; a[i-1] = tmp;
}
}
static void pt_solve1_perturb(const pt_svopt_t *opt, const pt_match_t *ma, uint64_t *x, solve_aux_t *aux)
{
uint32_t i;
double y;
for (i = 0; i < aux->size; ++i) {
uint32_t k = (uint32_t)ma->cc[aux->off + i];
y = kr_drand_r(x);
if (y < opt->f_perturb)
aux->s[k] = -aux->s[k];
}
}
static int64_t pt_solve1_optimize(const pt_match_t *ma, uint32_t topn, uint64_t *x, solve_aux_t *aux, uint32_t *n_iter)
{
uint32_t i;
while (1) {
uint32_t n_flip = 0;
++(*n_iter);
ks_shuffle_uint32_t(aux->size, aux->shuffled, x);
for (i = 0; i < aux->size; ++i) {
uint32_t k = aux->shuffled[i];
uint32_t o = ma->idx.a[k] >> 32;
uint32_t n = (uint32_t)ma->idx.a[k], j;
uint64_t z[2];
int8_t s;
for (j = 0; j < n; ++j) {
const pt_match1_t *m = &ma->ma.a[o + j];
///assert(m->sid[0] == k);
aux->buf[j] = (uint64_t)((uint32_t)-1 - m->w) << 32 | (o + j);
}
radix_sort_gfa64(aux->buf, aux->buf + n);///still sort by edge weight
for (j = 0, z[0] = z[1] = 0; j < n && j < topn; ++j) {
const pt_match1_t *m = &ma->ma.a[(uint32_t)aux->buf[j]];
if (aux->s[m->sid[1]] > 0) z[0] += m->w;
else if (aux->s[m->sid[1]] < 0) z[1] += m->w;
}
if (z[0] == z[1]) continue;
s = z[0] > z[1]? -1 : 1;
if (aux->s[k] != s)
aux->s[k] = s, ++n_flip;
}
if (n_flip == 0) break;
}
return pt_score(ma, topn, aux);
}
uint32_t pt_solve1(const pt_svopt_t *opt, const pt_match_t *ma, uint64_t *x, solve_aux_t *aux)
{
uint32_t j, k, n_iter = 0;
int64_t sc_ori, sc_opt = -(1<<30), sc;
if (aux->size < 2) return 0;///how many unitigs
// first guess
///randomly assign haplotype status, and get a score
sc_ori = pt_solve1_init_phase(ma, opt->topn, x, aux);
if (aux->size == 2) return 0;
// optimize
sc_opt = pt_solve1_optimize(ma, opt->topn, x, aux, &n_iter);
for (j = 0; j < aux->size; ++j)
aux->s_tmp[aux->shuffled[j]] = aux->s[aux->shuffled[j]];
for (k = 0; k < (uint32_t)opt->n_perturb; ++k) {
pt_solve1_perturb(opt, ma, x, aux);
sc = pt_solve1_optimize(ma, opt->topn, x, aux, &n_iter);
if (sc > sc_opt) {
for (j = 0; j < aux->size; ++j)
aux->s_tmp[aux->shuffled[j]] = aux->s[aux->shuffled[j]];
sc_opt = sc;
} else {
for (j = 0; j < aux->size; ++j)
aux->s[aux->shuffled[j]] = aux->s_tmp[aux->shuffled[j]];
}
}
for (j = 0; j < aux->size; ++j)
aux->s[aux->shuffled[j]] = aux->s_tmp[aux->shuffled[j]];
fprintf(stderr, "[%s] group:%d, size:%d, #edges:%u, #iter:%d, sc_ori:%ld, sc_opt:%ld\n", __func__,
(uint32_t)(ma->cc[aux->off]>>32), (uint32_t)(aux->size), (uint32_t)(aux->a.n), n_iter, (long)sc_ori, (long)sc_opt);
return n_iter;
}
int8_t *pt_solve_core(const pt_svopt_t *opt, const pt_match_t *ma)
{
int8_t *s;
uint32_t st, i, max = 0;
uint64_t x = opt->seed;
solve_aux_t *aux;
CALLOC(aux, 1); CALLOC(aux->s, ma->n_seq); CALLOC(aux->s_tmp, ma->n_seq);
kv_init(aux->a);
for (i = 0; i < ma->n_seq; ++i) {///count how many links for each unitig
uint32_t n = pt_n(*ma, i);
max = max > n? max : n;
}
MALLOC(aux->buf, max);
MALLOC(aux->shuffled, ma->n_seq); // FIXME: this is over-allocation for convenience
for (st = 0, i = 1; i <= ma->n_seq; ++i) {
if (i == ma->n_seq || ma->cc[st]>>32 != ma->cc[i]>>32) {///at same group
if (i - st >= 2) {///all unitigs in the same group
aux->off = st, aux->size = i - st;
pt_solve1(opt, ma, &x, aux);
}
st = i;
}
}
s = aux->s;
kv_destroy(aux->a); free(aux->buf); free(aux->shuffled); free(aux->s_tmp);
free(aux);
return s;
}
void set_trio_flag(ma_ug_t *ug, asg_t *read_g, uint32_t uID, uint8_t* trio_flag, trans_chain* t_ch,
hap_overlaps_list* ha, pt_match_t *ma, int8_t s)
{
uint32_t i;
ma_utg_t *u = &(ug->u.a[uID]);
for (i = 0; i < u->n; i++)
{
trio_flag[u->a[i]>>33] |= SET_TRIO;
if(t_ch->is_r_het[u->a[i]>>33] == N_HET) continue;
if(s == 0)
{
if(t_ch->is_r_het[u->a[i]>>33]&P_HET)//special case
{
trio_flag[u->a[i]>>33] |= FATHER;
}
continue;
}
trio_flag[u->a[i]>>33] |= (s > 0? FATHER:MOTHER);
}
// pt_match1_t *o = pt_a(*ma, uID);
// uint32_t n = pt_n(*ma, uID), k, qn, tn, qs, qe, r_beg, r_end, offset, oLen, found;
// int index;
// for (k = 0; k < n; ++k)
// {
// qn = o[k].sid[0]; tn = o[k].sid[1]; qs = qe = (uint32_t)-1;
// index = get_specific_hap_overlap(&(ha->x[qn]), qn, tn);
// if(index != -1 && ha->x[qn].a.a[index].score == (long long)o[k].w)
// {
// qs = ha->x[qn].a.a[index].x_beg_pos;
// qe = ha->x[qn].a.a[index].x_end_pos - 1;
// }
// else
// {
// index = get_specific_hap_overlap(&(ha->x[tn]), tn, qn);
// if(index != -1 && ha->x[tn].a.a[index].score == (long long)o[k].w)
// {
// qs = ha->x[qn].a.a[index].y_beg_pos;
// qe = ha->x[qn].a.a[index].y_end_pos - 1;
// }
// }
// if(qs == (uint32_t)-1 || qe == (uint32_t)-1) fprintf(stderr, "ERROR\n");
// for (i = 0, offset = 0, found = 0; i < u->n; i++)
// {
// r_beg = offset; r_end = offset + (long long)(read_g->seq[u->a[i]>>33].len) - 1;
// offset += (uint32_t)u->a[i];
// oLen = ((MIN(qe, r_end) >= MAX(qs, r_beg))? MIN(qe, r_end) - MAX(qs, r_beg) + 1 : 0);
// if(oLen > 0) found = 1;
// if(found == 1 && oLen == 0) break;
// if(oLen > 0 && t_ch->is_r_het[u->a[i]>>33] != N_HET)
// {
// trio_flag[u->a[i]>>33] |= (s > 0? FATHER:MOTHER);
// }
// }
// }
}
void pt_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag)
{
pt_svopt_t opt;
int8_t *s = NULL;
uint64_t *buf = NULL, i;
pt_svopt_init(&opt);
pt_g_t *pg = init_pt_g_t(ovlp, t_ch, ug, read_g, f_rate);
///debug_pt_g_t(pg);
pt_cc(pg->e);
s = pt_solve_core(&opt, pg->e);
MALLOC(buf, pg->e->ma.n); // FIXME: this is over-allocation for convenience
for (i = 0; i < pg->e->n_seq; ++i) {
uint64_t z[2];
uint32_t o = pg->e->idx.a[i] >> 32;
uint32_t n = (uint32_t)pg->e->idx.a[i], j;
set_trio_flag(ug, read_g, i, trio_flag, t_ch, ovlp, pg->e, s[i]);
pg->info.a[i].s = s[i];
for (j = 0; j < n; ++j) {
const pt_match1_t *m = &pg->e->ma.a[o + j];
buf[j] = (uint64_t)((uint32_t)-1 - m->w) << 32 | (o + j);
}
radix_sort_gfa64(buf, buf + n);
for (j = 0, z[0] = z[1] = 0; j < n; ++j) {
const pt_match1_t *m = &pg->e->ma.a[(uint32_t)buf[j]];
if (s[m->sid[1]] > 0) z[0] += m->w;
else if (s[m->sid[1]] < 0) z[1] += m->w;
}
pg->info.a[i].m[0] = z[0], pg->info.a[i].m[1] = z[1];
}
free(buf);
free(s);
destory_pt_g_t(&pg);
}
partig.h
+51
View File
@@ -0,0 +1,51 @@
#ifndef __PURTIG__
#define __PURTIG__
#include <stdio.h>
#include <stdint.h>
#include "kvec.h"
#include "Overlaps.h"
#include "Purge_Dups.h"
typedef struct {
uint32_t bS, bE;
uint32_t nS, nE;
uint32_t uID;
uint8_t hs;
}pt_node_t;
typedef struct {
///sid[0]: query id
///sid[1]: target id
uint32_t sid[2];
uint32_t w;
} pt_match1_t;
typedef struct {
kvec_t(uint64_t) idx;
kvec_t(pt_match1_t) ma;
uint64_t* cc;
uint32_t n_seq;
} pt_match_t;
typedef struct {
///cnt1: how many unique minimizers
///cnt2: how many non-unique minimizers
///uint32_t cnt2, cnt1;
uint64_t m[2];
int8_t s;
} pt_uinfo_t;
typedef struct {
kvec_t(pt_node_t) p;
ma_ug_t *ug;
asg_t *rg;
trans_chain* t_ch;
///kvec_t(int8_t) s; ///status
kvec_t(pt_uinfo_t) info; ///status
kvec_t(uint32_t) p_idx;
pt_match_t* e;
}pt_g_t;
void pt_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag);
#endif