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

code clean

Browse Source
This commit is contained in:
chhylp123
2022-11-17 09:31:56 -05:00
parent 2c7b164f11
commit 9831afaf26
12 changed files with 1296 additions and 127 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+1
View File
@@ -1775,6 +1775,7 @@ int ha_assemble(void)
{
// debug_mc_g_t(MC_NAME);
// debug_mc_gg_t(MC_NAME, 0, 0);
// quick_debug_phasing(MC_NAME);
extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
int r, hom_cov = -1, ovlp_loaded = 0;
if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
CommandLines.cpp
+6
View File
@@ -143,6 +143,12 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " --l-msjoin INT\n");
fprintf(stderr, " detect misjoined unitigs of >=INT in size; 0 to disable [%lu]\n", asm_opt->misjoin_len);
fprintf(stderr, " Ultra-Long-integration (beta):\n");
fprintf(stderr, " --ul FILEs file names of Ultra-Long reads [r1.fq,r2.fq,...]\n");
fprintf(stderr, " --ul-rate FLOAT\n");
fprintf(stderr, " similarity threshold for UL-to-HiFi alignment [%.3g]\n", asm_opt->ul_error_rate);
fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
fprintf(stderr, "See `https://hifiasm.readthedocs.io/en/latest/' or `man ./hifiasm.1' for complete documentation.\n");
}
CommandLines.h
+1 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.17.2-r433"
#define HA_VERSION "0.17.3-r439"
#define VERBOSE 0
Overlaps.cpp
+98 -50
View File
@@ -817,7 +817,7 @@ long long get_specific_overlap(ma_hit_t_alloc* x, uint32_t qn, uint32_t tn)
inline void set_reverse_overlap(ma_hit_t* dest, ma_hit_t* source)
void set_reverse_overlap(ma_hit_t* dest, ma_hit_t* source)
{
dest->qns = Get_tn(*source);
dest->qns = dest->qns << 32;
@@ -924,7 +924,7 @@ void normalize_ma_hit_t_single_side_advance(ma_hit_t_alloc* sources, long long n
}
}
// if(VERBOSE >= 1)
if(VERBOSE >= 1)
{
fprintf(stderr, "[M::%s] takes %0.2fs\n\n", __func__, Get_T()-startTime);
}
@@ -10140,7 +10140,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FIL
v = au[j].v;
fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\tL2:i:%u\n",
prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1],
prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), au[j].ou);
prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), 0/**au[j].ou**/);
}
@@ -10153,7 +10153,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FIL
v = au[j].v;
fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\tL2:i:%u\n",
prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1],
prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), au[j].ou);
prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), 0/**au[j].ou**/);
}
}
}
@@ -13674,14 +13674,14 @@ long long gap_fuzz, bub_label_t* b_mask_t)
if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(cov->t_ch, output_file_name);
}
///for debug
// char* gfa_name = (char*)malloc(strlen(output_file_name)+50);
// sprintf(gfa_name, "%s.pre.clean_d_utg.noseq.gfa", output_file_name);
// FILE* output_file = fopen(gfa_name, "w");
// ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file);
// fclose(output_file);
// free(gfa_name);
///for debug
hic_analysis(ug, sg, cov?cov->t_ch:t_ch, &opt, 0, asm_opt.scffold?&rhits:NULL);
@@ -16641,52 +16641,79 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
redo:
///print_untig((ug), 61955, "i-0:", 0);
// fprintf(stderr, "[M::%s] 0\n", __func__);
asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, NULL, 1);
// 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__);
/**********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__);
/**********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__);
pre_cons = get_graph_statistic(g);
// fprintf(stderr, "[M::%s] 5\n", __func__);
///need consider tangles
asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, NULL, 1);
// fprintf(stderr, "[M::%s] 6\n", __func__);
/**********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);
// 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);
// 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);
// 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);
// 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);
// fprintf(stderr, "[M::%s] 11\n", __func__);
///need consider tangles
///note we need both the read graph and the untig graph
}
/**********debug**********/
cur_cons = get_graph_statistic(g);
// fprintf(stderr, "[M::%s] 12\n", __func__);
}
if(just_bubble_pop == 0)
{
// fprintf(stderr, "[M::%s] 13\n", __func__);
cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
// 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__);
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__);
// 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);
// fprintf(stderr, "[M::%s] 17\n", __func__);
drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het);
// 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);
// 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);
// fprintf(stderr, "[M::%s] 20\n", __func__);
goto redo;
}
}
@@ -17624,12 +17651,13 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
__func__, asm_opt.recover_atg_cov_min);
}
// fprintf(stderr, "[M::%s] 0\n", __func__);
adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t, cov->is_r_het);
// fprintf(stderr, "[M::%s] 1\n", __func__);
///primary_flag = get_utg_attributes(*ug, read_g, coverage_cut, sources, ruIndex);
update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, cov->is_r_het, 0,
DOUBLE_CHECK_THRES, flag, drop_rate);
// fprintf(stderr, "[M::%s] 2\n", __func__);
nsg = (*ug)->g;
n_vtx = nsg->n_seq;
for (v = 0; v < n_vtx; ++v)
@@ -17638,36 +17666,44 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
nsg->seq[v].c = PRIMARY_LABLE;
EvaluateLen((*ug)->u, v) = (*ug)->u.a[v].n;
}
// fprintf(stderr, "[M::%s] 3\n", __func__);
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);
// fprintf(stderr, "[M::%s] 4\n", __func__);
///delete_useless_nodes(ug);
delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex);
// fprintf(stderr, "[M::%s] 5\n", __func__);
update_hap_label(*ug, read_g);
// fprintf(stderr, "[M::%s] 6\n", __func__);
update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, cov->is_r_het, 0,
DOUBLE_CHECK_THRES, flag, drop_rate);
// fprintf(stderr, "[M::%s] 7\n", __func__);
force_trio_clean((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, flag, 0.55, 0.01, 5);
// fprintf(stderr, "[M::%s] 8\n", __func__);
///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);
// fprintf(stderr, "[M::%s] 9\n", __func__);
if (!(asm_opt.flag & HA_F_BAN_POST_JOIN))
{
rescue_missing_overlaps_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang,
min_ovlp, 0, 1, NULL, b_mask_t);
// fprintf(stderr, "[M::%s] 10\n", __func__);
renew_utg(ug, read_g, new_rtg_edges);
// fprintf(stderr, "[M::%s] 11\n", __func__);
rescue_contained_reads_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang,
min_ovlp, 10, 0, 1, NULL, NULL, b_mask_t);
// fprintf(stderr, "[M::%s] 12\n", __func__);
renew_utg(ug, read_g, new_rtg_edges);
// fprintf(stderr, "[M::%s] 13\n", __func__);
}
///if(flag == MOTHER) print_untig_by_read(*ug, "m64043_200627_000137/124716590/ccs", 2789716, NULL, NULL, "beg");
@@ -17675,13 +17711,17 @@ 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, cov->is_r_het, 1,
FINAL_DOUBLE_CHECK_THRES, flag, drop_rate);
// fprintf(stderr, "[M::%s] 14\n", __func__);
update_hap_label(NULL, read_g);
// fprintf(stderr, "[M::%s] 15\n", __func__);
renew_utg(ug, read_g, new_rtg_edges);
// fprintf(stderr, "[M::%s] 16\n", __func__);
///delete_useless_nodes(ug);
delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex);
// fprintf(stderr, "[M::%s] 17\n", __func__);
if(asm_opt.purge_level_trio == 1)
@@ -17692,13 +17732,15 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
///delete_useless_nodes(ug);
delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex);
}
// fprintf(stderr, "[M::%s] 18\n", __func__);
set_drop_trio_flag(*ug);
// fprintf(stderr, "[M::%s] 19\n", __func__);
destory_hap_cov_t(&cov);
// fprintf(stderr, "[M::%s] 20\n", __func__);
// purge_dump(*ug);
renew_utg(ug, read_g, new_rtg_edges);
// fprintf(stderr, "[M::%s] 21\n", __func__);
}
@@ -17736,28 +17778,23 @@ char *f_prefix, uint8_t *kpt_buf, kvec_asg_arc_t_warp *r_edges)
kvec_asg_arc_t_warp new_rtg_edges;
kv_init(new_rtg_edges.a);
adjust_utg_by_trio(&ug, sg, flag, 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, b_mask_t);
if(asm_opt.b_low_cov > 0)
{
break_ug_contig(&ug, sg, &R_INF, coverage_cut, sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp,
&asm_opt.b_low_cov, NULL, asm_opt.m_rate);
}
if(asm_opt.b_high_cov > 0)
{
break_ug_contig(&ug, sg, &R_INF, coverage_cut, sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp,
NULL, &asm_opt.b_high_cov, asm_opt.m_rate);
}
if(kpt_buf)
{
update_dump_trio(R_INF.trio_flag, sg->n_seq, kpt_buf, ug);
}
if(is_bench)
{
free(gfa_name);
@@ -17775,7 +17812,6 @@ char *f_prefix, uint8_t *kpt_buf, kvec_asg_arc_t_warp *r_edges)
///debug_untig_length(ug, tipsLen, gfa_name);
///print_untig_by_read(ug, "m64011_190901_095311/125831121/ccs", 2310925, "end");
ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 1);
ma_ug_print(ug, sg, coverage_cut, sources, ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file);
fclose(output_file);
@@ -21075,7 +21111,7 @@ uint32_t type)
kvec_t(uint32_t) u_vecs;
uint32_t i = 0, maxEvaluateLen = 0, maxBaseLen = 0, v, totalEvaluateLen = 0;
u_vecs.a = a, u_vecs.n = n;
if(u_vecs.n > 0)
if(u_vecs.n > 0) ///u_vecs does not contain startID && endId
{
v_x = &(ug->u.a[u_vecs.a[0]>>1]);
asg_seq_del(nsg, u_vecs.a[0]>>1);
@@ -23476,7 +23512,6 @@ buf_t* bb, uint8_t* visit, uint32_t trio_flag, long long totalNodeLen, long long
}
n_reduce = drop_useless_edges(ug, bb, visit, beg, end, nodes, nodes_n, totalNodeLen);
}
// if(pop_bubble_at_tangle(ug, 10000000, nodes, nodes_n, beg, end, trio_flag, DROP)!=0)
// {
// fprintf(stderr, "false bubble popping: beg: %u, end: %u\n", beg>>1, end>>1);
@@ -23497,19 +23532,18 @@ buf_t* bb, uint8_t* visit, uint32_t trio_flag, long long totalNodeLen, long long
v = w;
if(v == beg) break;
}
return is_found;
}
uint32_t cut_edges_progressive(ma_ug_t *ug, kvec_t_u64_warp* edges, float drop_ratio,
uint64_t* nodes, uint64_t nodes_n, uint32_t beg, uint32_t end, buf_t* bb, uint8_t* visit,
uint32_t trio_flag, long long totalNodeLen, long long totalBaseLen)
uint32_t trio_flag, long long totalNodeLen, long long totalBaseLen, uint32_t max_arc_n)
{
uint32_t k, v, i, kv, w, nv, ov_max, ban, is_found = 0;
asg_arc_t *a = NULL, *av = NULL;
asg_t* nsg = ug->g;
radix_sort_arch64(edges->a.a, edges->a.a + edges->a.n);
for (k = 0; k < edges->a.n; k++)
for (k = 0; k < edges->a.n && k < max_arc_n; k++)
{
a = &nsg->arc[(uint32_t)edges->a.a[k]];
if(a->del) continue;
@@ -23586,20 +23620,17 @@ float drop_ratio)
// if(debug_is_circle==0) fprintf(stderr, "Not circle: beg: %u, end: %u\n", beg>>1, end>>1);
// if(debug_is_circle==1) fprintf(stderr, "Circle: beg: %u, end: %u\n", beg>>1, end>>1);
/**************************debug**************************/
is_found = process_tangles(ug, nodes, nodes_n, beg, end, bb, visit, trio_flag, totalNodeLen,
totalBaseLen);
///if(is_found == 1) fprintf(stderr, "***Found: beg>>1: %u, end>>1: %u\n", beg>>1, end>>1);
if(is_found == 0)
{
is_found = cut_edges_progressive(ug, edges, drop_ratio, nodes, nodes_n,
beg, end, bb, visit, trio_flag, totalNodeLen, totalBaseLen);
beg, end, bb, visit, trio_flag, totalNodeLen, totalBaseLen, 48);
///if(is_found == 1) fprintf(stderr, "***Cutting Found: beg>>1: %u, end>>1: %u\n", beg>>1, end>>1);
}
if(is_found == 1)
{
@@ -23631,10 +23662,8 @@ float drop_ratio)
is_found = 0;
}
}
recover_edges(nsg, edges, nodes, nodes_n, beg, end, beg_c, end_c, 1-is_found);
}
void resolve_tangles(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig,
@@ -23674,9 +23703,8 @@ uint32_t trio_flag, float drop_ratio)
uint8_t* visit = NULL;
visit = (uint8_t*)malloc(sizeof(uint8_t) * nsg->n_seq);
uint32_t n_reduce, flag;
uint32_t n_reduce, flag, dbg_round = 0;
n_vtx = nsg->n_seq * 2;
while (1)
{
n_reduce = 0;
@@ -23699,8 +23727,8 @@ uint32_t trio_flag, float drop_ratio)
}
}
if(n_reduce == 0) break;
dbg_round++;
}
asg_cleanup(nsg);
asg_symm(nsg);
@@ -23720,7 +23748,6 @@ uint32_t trio_flag, float drop_ratio)
}
}
n_vtx = nsg->n_seq;
for (v = 0; v < n_vtx; ++v)
{
@@ -23740,7 +23767,6 @@ uint32_t trio_flag, float drop_ratio)
}
nsu->n = m;
}
n_vtx = nsg->n_seq;
for (i = 0; i < n_vtx; ++i)
{
@@ -23764,7 +23790,6 @@ uint32_t trio_flag, float drop_ratio)
if(get_real_length(nsg, w^1, NULL) != 1) continue;
if(IsMerge(ug->u, w>>1) != 0) continue;
end = w;
///we have three types of merged nodes
///1) CONVEX_M: one direction has two out-nodes, another direction has one out-node
///2) UNROLL_E: one direction has one out-node, another direction doesn't has out-node
@@ -23774,13 +23799,11 @@ uint32_t trio_flag, float drop_ratio)
unroll_tangle(src, ug, nsu->a, nsu->n, beg, end, &e_vecs, trio_flag, &b_0, visit, drop_ratio);
}
for (v = 0; v < ruIndex->len; v++)
{
get_R_to_U(ruIndex, v, &uId, &is_Unitig);
if(is_Unitig == 1) ruIndex->index[v] = (uint32_t)-1;
}
kv_destroy(u_vecs.a);
kv_destroy(e_vecs.a);
@@ -23799,7 +23822,6 @@ uint32_t trio_flag, float drop_ratio)
EvaluateLen(src->u, v) = src->u.a[v].n;
}
///print_untig_by_read(src, "m64076_200203_181219/82511682/ccs", 2429597, NULL, NULL, "end-1");
}
@@ -24842,18 +24864,20 @@ int load_ruIndex(R_to_U* ruIndex, char* read_file_name)
int f_flag = 0;
f_flag += fread(&(ruIndex)->len, sizeof((ruIndex)->len), 1, fp);
(ruIndex)->index = (uint32_t*)malloc(sizeof(uint32_t)*(ruIndex)->len);
f_flag += fread((ruIndex)->index, sizeof((ruIndex)->index[0]), (ruIndex)->len, fp);
f_flag += fread((ruIndex)->index, sizeof((*((ruIndex)->index))), (ruIndex)->len, fp);
if(!(asm_opt.ar)) {
R_INF.trio_flag = (uint8_t*)malloc(sizeof(uint8_t)*(ruIndex)->len);
f_flag += fread(R_INF.trio_flag, sizeof(R_INF.trio_flag[0]), (ruIndex)->len, fp);
f_flag += fread(R_INF.trio_flag, sizeof((*(R_INF.trio_flag))), (ruIndex)->len, fp);
} else {
fseek(fp, sizeof((*(R_INF.trio_flag)))*(ruIndex)->len, SEEK_CUR);
}
// CALLOC(ruIndex->is_het, ruIndex->len);
// f_flag += fread(ruIndex->is_het, 1, ruIndex->len, fp);
// f_flag += fread(&(asm_opt.hom_global_coverage_set), sizeof(asm_opt.hom_global_coverage_set), 1, fp);
// f_flag += fread(&(asm_opt.hom_global_coverage), sizeof(asm_opt.hom_global_coverage), 1, fp);
f_flag += fread(&(asm_opt.hom_global_coverage_set), sizeof(asm_opt.hom_global_coverage_set), 1, fp);
f_flag += fread(&(asm_opt.hom_global_coverage), sizeof(asm_opt.hom_global_coverage), 1, fp);
free(index_name);
fflush(fp);
@@ -29291,9 +29315,8 @@ int max_hang, int min_ovlp, bubble_type* bub, long long gap_fuzz)
void rescue_bubble_by_chain(asg_t *sg, ma_sub_t *coverage_cut, 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, uint32_t chainLenThres, long long gap_fuzz,
bub_label_t* b_mask_t)
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, uint32_t chainLenThres, long long gap_fuzz, bub_label_t* b_mask_t, long long no_trio_recover)
{
kvec_asg_arc_t_warp new_rtg_edges;
kv_init(new_rtg_edges.a);
@@ -29327,7 +29350,7 @@ bub_label_t* b_mask_t)
beg_idx = bub.f_bub; occ = bub.b_bub + bub.b_end_bub + bub.tangle_bub;
rescue_bubbles_by_missing_ovlp_backward(ug, sg, sources, coverage_cut, ruIndex, max_hang, min_ovlp, chainLenThres, beg_idx, occ, &bub, b_mask_t);
if(ha_opt_triobin(&asm_opt))
if((!no_trio_recover) && (ha_opt_triobin(&asm_opt)))
{
ma_ug_destroy(ug); ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
reset_bub(&bub, ug, cov->t_ch, &new_rtg_edges);
@@ -31675,6 +31698,25 @@ ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t,
return sg;
}
void renew_g(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, long long *n_read,
uint64_t **readLen, ma_sub_t **coverage_cut, R_to_U *ruIndex, asg_t **sg,
int64_t mini_overlap_length, int64_t max_hang_length,
ug_opt_t *uopt, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio,
int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file)
{
ma_ug_t *iug = ul_realignment_gfa(uopt, *sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
asm_opt.max_short_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file);
asg_t *ng = gen_ng(iug, *sg, uopt, coverage_cut, ruIndex, 256);
ma_ug_destroy(iug); asg_destroy(*sg);
(*sources) = R_INF.paf;
(*reverse_sources) = R_INF.reverse_paf;
(*n_read) = R_INF.total_reads;
(*readLen) = R_INF.read_length;
(*sg) = ng;
ma_hit_contained_advance(*sources, *n_read, *coverage_cut, ruIndex, max_hang_length, mini_overlap_length);
post_rescue(uopt, *sg, (*sources), (*reverse_sources), ruIndex, b_mask_t, 0);
}
void clean_graph(
int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
long long n_read, uint64_t* readLen, long long mini_overlap_length,
@@ -31770,11 +31812,17 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
debug_gfa:;
gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex,
(asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, &b_mask_t);
set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length);
// set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length);
}
if(asm_opt.ar) {
ul_realignment_gfa(&uopt, sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
asm_opt.max_short_tip, &b_mask_t, ha_opt_triobin(&asm_opt), o_file);
renew_g(&sources, &reverse_sources, &n_read, &readLen, &coverage_cut, ruIndex, &sg, mini_overlap_length, max_hang_length,
&uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, asm_opt.max_short_tip, &b_mask_t,
ha_opt_triobin(&asm_opt), o_file);
// ma_ug_t *iug = ul_realignment_gfa(&uopt, sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
// asm_opt.max_short_tip, &b_mask_t, ha_opt_triobin(&asm_opt), o_file);
// gen_ng(iug, sg, &uopt, &coverage_cut, ruIndex, 100);
// exit(1);
}
// print_debug_gfa(sg, NULL, coverage_cut, "UL.debug", sources, ruIndex, max_hang_length, mini_overlap_length, 0, 0, 0);
/**
Overlaps.h
+4 -3
View File
@@ -32,7 +32,7 @@
// #define PRIMARY_LABLE 1
// #define ALTER_LABLE 2
// #define HAP_LABLE 4
#define HA_RE_UL_ID "re"
#define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32))
#define Get_qs(RECORD) ((uint32_t)((RECORD).qns))
@@ -880,7 +880,7 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, int m
void rescue_bubble_by_chain(asg_t *sg, ma_sub_t *coverage_cut, 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, uint32_t chainLenThres, long long gap_fuzz,
bub_label_t* b_mask_t);
bub_label_t* b_mask_t, long long no_trio_recover);
typedef struct{
double weight;
@@ -1140,7 +1140,8 @@ void hic_clean(asg_t* read_g);
int64_t count_edges_v_w(asg_t *g, uint32_t v, uint32_t w);
void renew_utg(ma_ug_t **ug, asg_t* read_g, kvec_asg_arc_t_warp* edge);
void merge_unitig_content(ma_utg_t* collection, ma_ug_t* ug, asg_t* read_g, kvec_asg_arc_t_warp* edge);
void reset_bub_label_t(bub_label_t* x, asg_t *g, uint64_t bub_dist, uint32_t check_cross);
void set_reverse_overlap(ma_hit_t* dest, ma_hit_t* source);
// void break_ug_contig(ma_ug_t **ug, asg_t *read_g, All_reads *RNF, ma_sub_t *coverage_cut,
// ma_hit_t_alloc* sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp,
// int* b_low_cov, int* b_high_cov, double m_rate);
Purge_Dups.cpp
+1 -1
View File
@@ -5623,7 +5623,7 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle
if(asm_opt.polyploidy <= 2)
{
mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1);
mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1, 0);
}
if(collect_p_trans && collect_p_trans_f == 1)
gfa_ut.cpp
+918 -43
View File
File diff suppressed because it is too large Load Diff
gfa_ut.h
+3 -1
View File
@@ -23,11 +23,13 @@ void asg_arc_cut_bub_links(asg_t *g, asg64_v *in, float len_rat, float sec_len_r
void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float ou_rat, uint32_t is_ou, bub_label_t *b_mask_t);
uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_rat, uint32_t is_ou);
uint32_t asg_cut_semi_circ(asg_t *g, uint32_t lim_len, uint32_t is_clean);
void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio,
ma_ug_t *ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio,
double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file);
void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t max_hang, int64_t min_ovlp, int64_t ul_occ);
void normalize_gou(asg_t *g);
void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd);
asg_t *gen_ng(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, uint64_t scaffold_len);
void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover);
// void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail);
#endif
hic.cpp
+12 -7
View File
@@ -16497,7 +16497,7 @@ void optimize_u_trans(kv_u_trans_t *ovlp, kvec_pe_hit* hits, ha_ug_index* idx)
}
}
kt_u_trans_t_idx(&k_trans, idx->ug->g->n_seq);
mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1);
mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1, 0);
for (i = m = 0; i < ovlp->n; i++){
x = &(ovlp->a[i]);
if(x->del) continue;
@@ -16524,7 +16524,7 @@ ha_ug_index* idx, uint64_t test_block_flip, uint64_t n_perturb)
__func__, test_block_flip, n_perturb);
renew_kv_u_trans(k_trans, link, &sl->hits, &(idx->t_ch->k_trans), idx, bub, s->s, NULL, 0);
mc_solve(NULL, NULL, k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag,
(bub->round_id == 0? 1 : 0), s->s, 1, NULL, test_block_flip?&(idx->t_ch->k_trans):0, 0);
(bub->round_id == 0? 1 : 0), s->s, 1, NULL, test_block_flip?&(idx->t_ch->k_trans):0, 0, 0);
}
@@ -16596,6 +16596,9 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
// label_unitigs_sm(s->s, NULL, idx->ug);
// goto skip_flipping;
// }
// print_debug_gfa(idx->read_g, idx->ug, opt->coverage_cut, "hic.phasing", opt->sources, opt->ruIndex,
// opt->max_hang, opt->min_ovlp, 0, 0, 0);
s = init_ps_t(11, idx->ug->g->n_seq);
// debug_round_test(s, 11, &bub, &k_trans, &link, &sl, idx, 1000, 10000);
for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++)
@@ -16605,8 +16608,11 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
// if(bub.round_id == 0) init_phase(idx, &k_trans, &bub, s);
// update_trans_g(idx, &k_trans, &bub);
/*******************************for debug************************************/
// mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag,
// (bub.round_id == 0? 1 : 0), s->s, 1, (asm_opt.ar)?(&bub):(NULL), &(idx->t_ch->k_trans), 0,
// (((bub.round_id+1) == bub.n_round)?1:0));
mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag,
(bub.round_id == 0? 1 : 0), s->s, 1, /**&bub**/NULL, &(idx->t_ch->k_trans), 0);
(bub.round_id == 0? 1 : 0), s->s, 1, NULL, &(idx->t_ch->k_trans), 0, 0);
/*******************************for debug************************************/
label_unitigs_sm(s->s, NULL, idx->ug);
@@ -16639,9 +16645,8 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
///print_hc_links(&link, 0, &hap);
// print_kv_u_trans(&k_trans, &link, s->s);
///print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, idx->link, idx);
///print_hits(idx, &sl.hits, fn1);
// print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, NULL/**idx->link**/, idx);
// print_hits(idx, &sl.hits, fn1, fn2);
///print_debug_bubble_graph(&bub, idx->ug, asm_opt.output_file_name);
@@ -17369,4 +17374,4 @@ void hic_benchmark(ma_ug_t *ug, asg_t* read_g)
hic_short_align_bench(asm_opt.hic_reads[0], asm_opt.hic_reads[1], output_file_name, ug_index);
free(output_file_name);
}
}
inter.cpp
+21 -15
View File
@@ -4876,7 +4876,8 @@ void update_ul_vec_t_ug(const ul_idx_t *uref, ul_vec_t *rch, vec_mg_lchain_t *uc
rch->bb.a[rch->bb.a[k].pidx].aidx = k;
}
// fprintf(stderr, "+ulid->%ld\n", ulid);
uint32_t sp = (uint32_t)-1, ep = (uint32_t)-1; k = l/**a_n - 1**/;///start from the max chain
uint32_t sp = (uint32_t)-1, ep = (uint32_t)-1, ch_n = 0;
k = l/**a_n - 1**/;///start from the max chain
for (l = 0; k >= 0; ) {
if(sp == (uint32_t)-1 || rch->bb.a[k].qe <= sp) {
if(sp != (uint32_t)-1) l += ep - sp;
@@ -4886,6 +4887,7 @@ void update_ul_vec_t_ug(const ul_idx_t *uref, ul_vec_t *rch, vec_mg_lchain_t *uc
}
if(rch->bb.a[k].pidx == (uint32_t)-1) k = -1;
else k = rch->bb.a[k].pidx;
ch_n++;
}
rch->dd = 0;
@@ -4897,8 +4899,9 @@ void update_ul_vec_t_ug(const ul_idx_t *uref, ul_vec_t *rch, vec_mg_lchain_t *uc
rch->dd = 1;
} else if(l < ((int64_t)rch->rlen)*0.001) {
rch->dd = 2;
} else if(ch_n < rch->bb.n) {///multiple chain, might be useful for the scaffolding
rch->dd = 3;
}
// fprintf(stderr, "[M::%s::] rch->dd::%u, rch->bb.n::%u\n",
// __func__, rch->dd, (uint32_t)rch->bb.n);
}
@@ -8975,9 +8978,12 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call
// b->num_correct_base += b->correct.corrected_base;
// b->num_recorrect_base += b->round2.dumy.corrected_base;
if(UL_INF.a[s->id+i].dd) {
free(s->seq[i]); s->seq[i] = NULL; b->num_correct_base++;
if(UL_INF.a[s->id+i].dd == 1 || UL_INF.a[s->id+i].dd == 2) {
b->num_correct_base++;
}
if(UL_INF.a[s->id+i].dd != 3) {
free(s->seq[i]); s->seq[i] = NULL;
}
s->hab[tid]->num_read_base++;
// fprintf(stderr, "[M::%s] rid:%ld, dd:%u\n", __func__, s->id+i, UL_INF.a[s->id+i].dd);
// int64_t mem[6], mem_hab[6];
@@ -9364,9 +9370,9 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac
// s->hab[i] = ha_ovec_buf_init(NULL, 0, 0, 1);
s->hab[i] = ha_ovec_init(0, 0, 1);
}
fprintf(stderr, "[M::%s::Start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
// fprintf(stderr, "[M::%s::Start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
kt_for(p->n_thread, worker_for_ul_scall_alignment, s, s->n);
fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
// fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
///debug
/**
uint64_t i;
@@ -9408,7 +9414,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac
p->num_bases += s->num_bases;
p->num_corrected_bases += s->num_corrected_bases;
p->num_recorrected_bases += s->num_recorrected_bases;
fprintf(stderr, "[M::%s::dump_start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
// fprintf(stderr, "[M::%s::dump_start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
for (i = 0; i < p->n_thread; ++i) {
push_uc_block_t(s->uopt, &(s->ll[i].tk), s->seq, s->len, s->id);
free(s->ll[i].tk.a);
@@ -9420,7 +9426,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac
}
free(s->seq[i]);
}
fprintf(stderr, "[M::%s::dump_done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
// fprintf(stderr, "[M::%s::dump_done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
/**
for (i = 0; i < (uint64_t)s->n; ++i) {
///debug
@@ -9495,10 +9501,10 @@ static void *worker_ul_rescall_pipeline(void *data, int step, void *in) // callb
for (i = 0; i < p->n_thread; ++i) {
s->hab[i] = ha_ovec_init(0, 0, 1); s->buf[i] = mg_tbuf_init();
}
fprintf(stderr, "[M::%s::Start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
// fprintf(stderr, "[M::%s::Start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
kt_for(p->n_thread, worker_for_ul_rescall_alignment, s, s->n);
fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
get_utepdat_t_mem(s, 1);
// fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
// get_utepdat_t_mem(s, 1);
for (i = 0; i < p->n_thread; ++i) {
p->num_bases += s->hab[i]->num_read_base;
@@ -9514,10 +9520,10 @@ static void *worker_ul_rescall_pipeline(void *data, int step, void *in) // callb
return s;
} else if (step == 2) { // step 3: dump
utepdat_t *s = (utepdat_t*)in; int64_t i, rid;
fprintf(stderr, "[M::%s::dump_start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
// fprintf(stderr, "[M::%s::dump_start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
for (i = 0; i < s->n; ++i) {
rid = s->id + i;
if(UL_INF.a[rid].dd == 0 && p->ucr_s && p->ucr_s->flag == 1) {
if(UL_INF.a[rid].dd == 3 && p->ucr_s && p->ucr_s->flag == 1) {///for the scaffolding
assert(s->seq[i]);
// if(s->seq[i] == NULL) fprintf(stderr, "[M::%s::]rid->%ld, len->%lu\n", __func__, rid, s->len[i]);
///for debug interval
@@ -9526,7 +9532,7 @@ static void *worker_ul_rescall_pipeline(void *data, int step, void *in) // callb
// if(UL_INF.a[rid].dd) fprintf(stderr, "rid->%ld\n", rid);
free(s->seq[i]);
}
fprintf(stderr, "[M::%s::dump_done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
// fprintf(stderr, "[M::%s::dump_done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
free(s->len); free(s->seq); free(s);
}
return 0;
@@ -14752,7 +14758,7 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_c
mg_idxopt_t opt; uldat_t sl;
int32_t cutoff;
char* gfa_name = NULL; MALLOC(gfa_name, strlen(asm_opt.output_file_name)+50);
sprintf(gfa_name, "%s.re", asm_opt.output_file_name);
sprintf(gfa_name, "%s.%s", asm_opt.output_file_name, HA_RE_UL_ID);
init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov);
cutoff = REA_ALIGN_CUTOFF;
rcut.cpp
+229 -5
View File
@@ -634,7 +634,7 @@ mb_g_t *init_mb_g_t(mc_match_t* e, kv_u_trans_t *ref, uint32_t is_sys)
for (m = 0; m < n; m++)
{
tn = ma_y(o[m]);
tb = p->u->idx.a[tn]>>1;
tb = p->u->idx.a[tn]>>1;///tn is the unitig id; tb is the block id
if(tb == qb)
{
continue;
@@ -1901,7 +1901,7 @@ uint8_t *lock, bits_p *vis, uint64_t id, mc_bp_res* r)
reset_mc_bp_iter(bp, &i, id);
memset(vis->a, 0, vis->n);
max_bid = max_uid = (uint32_t)-1;
f_bid = i.bid; f_uid = i.uid;
f_bid = i.bid; f_uid = i.uid;///f_bid::bubble id, f_uid::unitig id
while (1)
{
@@ -2902,14 +2902,20 @@ void print_hap_s(int8_t *s, uint32_t sn)
}
}
void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov)
void dump_debug_phasing(const char* fn, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate,
uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref);
void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov, int is_dump)
{
if(is_dump) {
dump_debug_phasing(MC_NAME, ta, ug, read_g, f_rate, renew_s, s, is_sys, bub, ref);
}
mc_opt_t opt;
mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed);
mc_g_t *mg = init_mc_g_t(ug, read_g, s, renew_s);
update_mc_edges(mg, ovlp, ta, t_ch, f_rate, is_sys);
// fprintf(stderr, "[M::%s:: # edges: %u]\n", __func__, (uint32_t)mg->e->ma.n);
fprintf(stderr, "[M::%s:: # edges: %u]\n", __func__, (uint32_t)mg->e->ma.n);
mb_solve_core(&opt, mg, ref, is_sys);
///debug_mc_g_t(mg);
@@ -3788,4 +3794,222 @@ void mc_solve_general(kv_u_trans_t *ta, uint32_t un, kv_gg_status *s, uint16_t h
destory_mc_gg_t(&mg);
if(write_dump) write_mc_gg_dump(ta, un, s, hapN, MC_NAME);
}
}
void dump_kv_u_trans_t(kv_u_trans_t *z, FILE *fp)
{
fwrite(&(z->n), sizeof(z->n), 1, fp);
fwrite(z->a, sizeof((*(z->a))), z->n, fp);
fwrite(&z->idx.n, sizeof(z->idx.n), 1, fp);
fwrite(z->idx.a, sizeof((*(z->idx.a))), z->idx.n, fp);
}
void dump_asg_t(asg_t *g, FILE *fp)
{
uint32_t tmp, Len;
tmp = g->n_arc;
fwrite(&tmp, sizeof(tmp), 1, fp);
tmp = g->is_srt;
fwrite(&tmp, sizeof(tmp), 1, fp);
tmp = g->n_seq;
fwrite(&tmp, sizeof(tmp), 1, fp);
tmp = g->is_symm;
fwrite(&tmp, sizeof(tmp), 1, fp);
tmp = g->r_seq;
fwrite(&tmp, sizeof(tmp), 1, fp);
// Len = g->n_seq*2;
// fwrite(g->seq_vis, sizeof((*(g->seq_vis))), Len, fp);
Len = g->n_seq*2;
fwrite(g->idx, sizeof((*(g->idx))), Len, fp);
fwrite(g->arc, sizeof((*(g->arc))), g->n_arc, fp);
fwrite(g->seq, sizeof((*(g->seq))), g->n_seq, fp);
}
void dump_ma_ug_t(ma_ug_t *ug, FILE *fp)
{
ma_utg_t *u = NULL; uint32_t t, i;
fwrite(&(ug->u.n), sizeof(ug->u.n), 1, fp);
for (i = 0; i < ug->u.n; i++) {
u = &(ug->u.a[i]);
t = u->len;
fwrite(&t, sizeof(t), 1, fp);
t = u->circ;
fwrite(&t, sizeof(t), 1, fp);
fwrite(&(u->start), sizeof(u->start), 1, fp);
fwrite(&(u->end), sizeof(u->end), 1, fp);
fwrite(&(u->n), sizeof(u->n), 1, fp);
fwrite(u->a, sizeof(uint64_t), u->n, fp);
}
dump_asg_t(ug->g, fp);
}
void dump_bubble_type(bubble_type* bub, FILE *fp)
{
fwrite(&(bub->chain_weight.n), sizeof(bub->chain_weight.n), 1, fp);
fwrite(bub->chain_weight.a, sizeof((*(bub->chain_weight.a))), bub->chain_weight.n, fp);
fwrite(&(bub->list.n), sizeof(bub->list.n), 1, fp);
fwrite(bub->list.a, sizeof((*(bub->list.a))), bub->list.n, fp);
fwrite(&(bub->num.n), sizeof(bub->num.n), 1, fp);
fwrite(bub->num.a, sizeof((*(bub->num.a))), bub->num.n, fp);
fwrite(&(bub->pathLen.n), sizeof(bub->pathLen.n), 1, fp);
fwrite(bub->pathLen.a, sizeof((*(bub->pathLen.a))), bub->pathLen.n, fp);
dump_ma_ug_t(bub->b_ug, fp);
dump_asg_t(bub->b_g, fp);
dump_ma_ug_t(bub->ug, fp);
}
void dump_debug_phasing(const char* fn, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate,
uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref)
{
fprintf(stderr, "\n[M::%s]\n", __func__);
char *buf = (char*)calloc(strlen(fn) + 50, 1);
sprintf(buf, "%s.hic.dbg.dump.bin", fn);
FILE *fp = fopen(buf, "w");
dump_kv_u_trans_t(ta, fp);///ta
dump_ma_ug_t(ug, fp);///ug
dump_asg_t(read_g, fp);///read_g
fwrite(&f_rate, sizeof(f_rate), 1, fp);///f_rate
fwrite(&renew_s, sizeof(renew_s), 1, fp);///renew_s
fwrite(s, sizeof((*s)), ug->g->n_seq, fp);///s
fwrite(&is_sys, sizeof(is_sys), 1, fp);///is_sys
dump_bubble_type(bub, fp);///bub
dump_kv_u_trans_t(ref, fp);///ta
fclose(fp);
free(buf);
}
void load_kv_u_trans_t(kv_u_trans_t *z, FILE *fp)
{
fread(&(z->n), sizeof(z->n), 1, fp);
z->m = z->n; MALLOC(z->a, z->n);
fread(z->a, sizeof((*(z->a))), z->n, fp);
fread(&z->idx.n, sizeof(z->idx.n), 1, fp);
z->idx.m = z->idx.n; MALLOC(z->idx.a, z->idx.n);
fread(z->idx.a, sizeof((*(z->idx.a))), z->idx.n, fp);
}
void load_asg_t(asg_t *g, FILE *fp)
{
uint32_t tmp, Len;
fread(&tmp, sizeof(tmp), 1, fp);
g->n_arc = g->m_arc = tmp;
fread(&tmp, sizeof(tmp), 1, fp);
g->is_srt = tmp;
fread(&tmp, sizeof(tmp), 1, fp);
g->n_seq = g->m_seq = tmp;
fread(&tmp, sizeof(tmp), 1, fp);
g->is_symm = tmp;
fread(&tmp, sizeof(tmp), 1, fp);
g->r_seq = tmp;
// Len = g->n_seq*2;
// fwrite(g->seq_vis, sizeof((*(g->seq_vis))), Len, fp);
Len = g->n_seq*2;
CALLOC(g->idx, Len); fread(g->idx, sizeof((*(g->idx))), Len, fp);
CALLOC(g->arc, g->n_arc); fread(g->arc, sizeof((*(g->arc))), g->n_arc, fp);
CALLOC(g->seq, g->n_seq); fread(g->seq, sizeof((*(g->seq))), g->n_seq, fp);
}
void load_ma_ug_t(ma_ug_t *ug, FILE *fp)
{
ma_utg_t *u = NULL; uint32_t t, i;
fread(&(ug->u.n), sizeof(ug->u.n), 1, fp);
ug->u.m = ug->u.n; CALLOC(ug->u.a, ug->u.n);
for (i = 0; i < ug->u.n; i++) {
u = &(ug->u.a[i]);
fread(&t, sizeof(t), 1, fp); u->len = t;
fread(&t, sizeof(t), 1, fp); u->circ = t;
fread(&(u->start), sizeof(u->start), 1, fp);
fread(&(u->end), sizeof(u->end), 1, fp);
fread(&(u->n), sizeof(u->n), 1, fp);
u->m = u->n; MALLOC(u->a, u->n);
fread(u->a, sizeof((*(u->a))), u->n, fp);
}
CALLOC(ug->g, 1);
load_asg_t(ug->g, fp);
}
void load_bubble_type(bubble_type* bub, FILE *fp)
{
fread(&(bub->chain_weight.n), sizeof(bub->chain_weight.n), 1, fp);
bub->chain_weight.m = bub->chain_weight.n; CALLOC(bub->chain_weight.a, bub->chain_weight.n);
fread(bub->chain_weight.a, sizeof((*(bub->chain_weight.a))), bub->chain_weight.n, fp);
fread(&(bub->list.n), sizeof(bub->list.n), 1, fp);
bub->list.m = bub->list.n; CALLOC(bub->list.a, bub->list.n);
fread(bub->list.a, sizeof((*(bub->list.a))), bub->list.n, fp);
fread(&(bub->num.n), sizeof(bub->num.n), 1, fp);
bub->num.m = bub->num.n; CALLOC(bub->num.a, bub->num.n);
fread(bub->num.a, sizeof((*(bub->num.a))), bub->num.n, fp);
fread(&(bub->pathLen.n), sizeof(bub->pathLen.n), 1, fp);
bub->pathLen.m = bub->pathLen.n; CALLOC(bub->pathLen.a, bub->pathLen.n);
fread(bub->pathLen.a, sizeof((*(bub->pathLen.a))), bub->pathLen.n, fp);
CALLOC(bub->b_ug, 1); load_ma_ug_t(bub->b_ug, fp);
CALLOC(bub->b_g, 1); load_asg_t(bub->b_g, fp);
CALLOC(bub->ug, 1); load_ma_ug_t(bub->ug, fp);
}
void load_debug_phasing(const char* fn, kv_u_trans_t **ta, ma_ug_t **ug, asg_t **read_g, double *f_rate,
uint32_t *renew_s, int8_t **s, uint32_t *is_sys, bubble_type **bub, kv_u_trans_t **ref)
{
fprintf(stderr, "\n[M::%s]\n", __func__);
char *buf = (char*)calloc(strlen(fn) + 50, 1);
sprintf(buf, "%s.hic.dbg.dump.bin", fn);
FILE *fp = fopen(buf, "r");
CALLOC((*ta), 1); load_kv_u_trans_t(*ta, fp);///ta
CALLOC((*ug), 1); load_ma_ug_t(*ug, fp);///ug
CALLOC((*read_g), 1); load_asg_t((*read_g), fp);///read_g
fread(f_rate, sizeof((*f_rate)), 1, fp);///f_rate
fread(renew_s, sizeof((*renew_s)), 1, fp);///renew_s
CALLOC((*s), (*ug)->g->n_seq); fread(*s, sizeof((*(*s))), (*ug)->g->n_seq, fp);///s
fread(is_sys, sizeof((*is_sys)), 1, fp);///is_sys
CALLOC((*bub), 1); load_bubble_type(*bub, fp);///bub
CALLOC((*ref), 1); load_kv_u_trans_t(*ref, fp);///ta
fclose(fp);
free(buf);
}
void quick_debug_phasing(const char* fn)
{
kv_u_trans_t *ta; ma_ug_t *ug; asg_t *read_g; double f_rate;
uint32_t renew_s; int8_t *s; uint32_t is_sys, k; bubble_type *bub; kv_u_trans_t *ref;
load_debug_phasing(fn, &ta, &ug, &read_g, &f_rate, &renew_s, &s, &is_sys, &bub, &ref);
mc_solve(NULL, NULL, ta, ug, read_g, f_rate, NULL, renew_s, s, is_sys, bub, ref, 0, 0);
for (k = 0; k < ug->g->n_seq; k++) {
fprintf(stderr, "utg%.6dl(len::%u)\n", (int32_t)(k)+1, ug->g->seq[k].len);
}
exit(1);
}
rcut.h
+2 -1
View File
@@ -118,11 +118,12 @@ static inline double kr_drand_r(uint64_t *x)
return u.d - 1.0;
}
void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov);
void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov, int is_dump);
void debug_mc_g_t(const char* name);
void mc_solve_general(kv_u_trans_t *ta, uint32_t un, kv_gg_status *s, uint16_t hapN, uint16_t update_ta, uint16_t write_dump);
kv_gg_status *init_mc_gg_status(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, R_to_U* ruIndex, uint64_t t_cov, uint16_t hapN);
void destory_mc_gg_t(mc_gg_t **p);
void debug_mc_gg_t(const char* fn, uint32_t update_ta, uint32_t convert_mc_g_t);
void quick_debug_phasing(const char* fn);
#endif