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

backup hic

Browse Source
This commit is contained in:
chhylp123
2021-04-18 21:57:16 -04:00
parent 49ead1ef35
commit 39c19618e8
2 changed files with 81 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Overlaps.cpp
+36 -27
View File
@@ -12578,7 +12578,7 @@ trans_chain* load_hc_hits(const char *fn)
fclose(fp);
free(buf);
// fprintf(stderr, "[M::%s::] ==> Hi-C cov have been loaded\n", __func__);
fprintf(stderr, "[M::%s::] ==> Hi-C cov have been loaded\n", __func__);
return t_ch;
}
@@ -12634,14 +12634,14 @@ bub_label_t* b_mask_t)
hic_analysis(ug, sg, cov?cov->t_ch:t_ch);
/**
char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
sprintf(gfa_name, "%s.d_utg.noseq.gfa", output_file_name);
FILE* output_file = fopen(gfa_name, "w");
ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file);
fclose(output_file);
free(gfa_name);
**/
// char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
// sprintf(gfa_name, "%s.d_utg.noseq.gfa", output_file_name);
// FILE* output_file = fopen(gfa_name, "w");
// ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file);
// fclose(output_file);
// free(gfa_name);
if(cov) destory_hap_cov_t(&cov);
if(t_ch) destory_trans_chain(&t_ch);
@@ -16210,6 +16210,8 @@ void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sou
fprintf(stderr, "bin files have been written.\n");
}
int load_debug_graph(asg_t** sg, ma_hit_t_alloc** sources, ma_sub_t** coverage_cut,
char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex);
int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, char* output_file_name)
{
char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
@@ -16218,6 +16220,14 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s
free(gfa_name);
return 0;
}
if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_debug_graph(NULL, NULL, NULL, output_file_name, NULL, NULL))
{
(*sources) = NULL;
(*reverse_sources) = NULL;
free(gfa_name);
return 1;
}
sprintf(gfa_name, "%s.ovlp.source", output_file_name);
if (!load_ma_hit_ts(sources, gfa_name)) {
free(gfa_name);
@@ -22858,8 +22868,7 @@ int write_ruIndex(R_to_U* ruIndex, char* read_file_name)
fwrite(&ruIndex->len, sizeof(ruIndex->len), 1, fp);
fwrite(ruIndex->index, sizeof(ruIndex->index[0]), ruIndex->len, fp);
fwrite(R_INF.trio_flag, sizeof(R_INF.trio_flag[0]), ruIndex->len, fp);
fwrite(ruIndex->is_het, 1, ruIndex->len, fp);
// fwrite(ruIndex->is_het, 1, ruIndex->len, fp);
free(index_name);
fflush(fp);
fclose(fp);
@@ -22884,8 +22893,8 @@ int load_ruIndex(R_to_U* ruIndex, char* read_file_name)
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);
CALLOC(ruIndex->is_het, ruIndex->len);
f_flag += fread(ruIndex->is_het, 1, ruIndex->len, fp);
// CALLOC(ruIndex->is_het, ruIndex->len);
// f_flag += fread(ruIndex->is_het, 1, ruIndex->len, fp);
free(index_name);
fflush(fp);
@@ -23099,7 +23108,11 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex)
fp = fopen(gfa_name, "r"); if(!fp) return 0;
sprintf(gfa_name, "%s.all.debug.asg_t.bin", output_file_name);
fp = fopen(gfa_name, "r"); if(!fp) return 0;
if((sg == NULL) || (sources == NULL) || (coverage_cut == NULL) || (reverse_sources == NULL) ||
(ruIndex == NULL))
{
return 1;
}
if((*sources)!=NULL)
{
@@ -29791,7 +29804,7 @@ float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name,
long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr,
ma_sub_t **coverage_cut_ptr, int debug_g)
{
char *o_file = NULL;
char *o_file = get_outfile_name(output_file_name);
ma_sub_t *coverage_cut = *coverage_cut_ptr;
asg_t *sg = *sg_ptr;
bub_label_t b_mask_t;
@@ -29982,24 +29995,20 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
ruIndex->is_het = rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3,
ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz, &b_mask_t);
if (asm_opt.flag & HA_F_VERBOSE_GFA)
{
/*******************************for debug***************************************/
write_debug_graph(sg, sources, coverage_cut, output_file_name, n_read, reverse_sources, ruIndex);
debug_gfa:;
/*******************************for debug***************************************/
set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex,
max_hang_length, mini_overlap_length);
}
o_file = get_outfile_name(output_file_name);
output_unitig_graph(sg, coverage_cut, o_file, sources, ruIndex, max_hang_length, mini_overlap_length);
flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL;
output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources,
asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length);
if (asm_opt.flag & HA_F_VERBOSE_GFA)
{
write_debug_graph(sg, sources, coverage_cut, output_file_name, n_read, reverse_sources, ruIndex);
debug_gfa:;
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))
{
if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION;
hic.cpp
+45 -6
View File
@@ -13514,6 +13514,39 @@ void idx_hc_links(kvec_pe_hit* hits, ha_ug_index* idx, bubble_type* bub)
}
}
inline uint32_t trans_checking_pass(bubble_type* bub, kv_u_trans_t *ref, uint32_t x, uint32_t y)
{
if(u_trans_n(*ref, x) == 0 || u_trans_n(*ref, y) == 0) return 0;
u_trans_t *a = NULL;
uint32_t n, i, f[2], qn, tn;
qn = x; tn = y;
a = u_trans_a(*ref, qn); n = u_trans_n(*ref, qn);
for (i = 0, f[0] = f[1] = 0; i < n; i++)
{
if(a[i].del) continue;
if(a[i].f == RC_2) continue;
if(IF_HOM(a[i].tn, *bub)) continue;
f[(a[i].tn == tn && a[i].f != RC_2)]++;
}
if(f[0] != 0) return 0;
if(f[1] == 0) return 0;
qn = y; tn = x;
a = u_trans_a(*ref, qn); n = u_trans_n(*ref, qn);
for (i = 0, f[0] = f[1] = 0; i < n; i++)
{
if(a[i].del) continue;
if(IF_HOM(a[i].tn, *bub)) continue;
f[(a[i].tn == tn && a[i].f != RC_2)]++;
}
if(f[0] != 0) return 0;
if(f[1] == 0) return 0;
fprintf(stderr, "M::%s::s-utg%.6ul<----->d-utg%.6ul\n", __func__, x+1, y+1);
return 1;
}
void weight_kv_u_trans(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub,
kv_u_trans_t *ta, trans_idx* dis)
{
@@ -13653,7 +13686,7 @@ pe_hit *hits, uint32_t occ, uint32_t qid, uint32_t qs, uint32_t qe, uint32_t tid
return w;
}
double get_hits_weight(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, trans_idx* dis,
double get_hits_weight(ha_ug_index* idx, bubble_type* bub, kvec_pe_hit* hits, hc_links* link, trans_idx* dis,
u_trans_t *t_a, uint32_t t_n, uint32_t qid, kv_u_trans_t *ta_idx)
{
/****************************may have bugs********************************/
@@ -13666,6 +13699,8 @@ u_trans_t *t_a, uint32_t t_n, uint32_t qid, kv_u_trans_t *ta_idx)
double w, i_w;
for (k = 0, w = 0; k < t_n; k++)
{
if(IF_HOM(qid, *bub)) continue;
if(IF_HOM(t_a[k].tn, *bub)) continue;
/****************************may have bugs********************************/
if(qid != t_a[k].tn)
{
@@ -13705,6 +13740,7 @@ u_trans_t *t_a, uint32_t t_n, uint32_t qid, kv_u_trans_t *ta_idx)
return w;
}
void adjust_weight_kv_u_trans(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub,
kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis)
{
@@ -13723,12 +13759,15 @@ kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis)
{
ta->a[m].nw /= (double)(MIN(hits->occ.a[ta->a[m].qn], hits->occ.a[ta->a[m].tn]));
}
/*******************************for debug************************************/
trans_checking_pass(bub, ref, ta->a[i].qn, ta->a[i].tn);
// if(trans_checking_pass(bub, ref, ta->a[i].qn, ta->a[i].tn)) ta->a[m].nw = 0;
/*******************************for debug************************************/
m++;
}
ta->n = m;
// fprintf(stderr, "+++++ta->n=%u\n", (uint32_t)ta->n);
for (k = 0; k < ta->idx.n; k++)///all nodes
{
if(IF_HOM(k, *bub)) continue;
@@ -13743,10 +13782,10 @@ kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis)
if(IF_HOM(a[i].qn, *bub)) continue;
if(IF_HOM(a[i].tn, *bub)) continue;
///(qn, tn^)
a[i].nw += get_hits_weight(idx, hits, link, dis,
a[i].nw += get_hits_weight(idx, bub, hits, link, dis,
u_trans_a(*ref, a[i].tn), u_trans_n(*ref, a[i].tn), a[i].qn, ta);
///(qn^1, tn)
a[i].nw += get_hits_weight(idx, hits, link, dis,
a[i].nw += get_hits_weight(idx, bub, hits, link, dis,
u_trans_a(*ref, a[i].qn), u_trans_n(*ref, a[i].qn), a[i].tn, ta);
}
@@ -13759,10 +13798,10 @@ kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis)
if(vis[tn]) continue;
if(qn == tn) continue;
///(qn, tn^)
w += get_hits_weight(idx, hits, link, dis,
w += get_hits_weight(idx, bub, hits, link, dis,
u_trans_a(*ref, tn), u_trans_n(*ref, tn), qn, ta);
///(qn^1, tn)
w += get_hits_weight(idx, hits, link, dis,
w += get_hits_weight(idx, bub, hits, link, dis,
u_trans_a(*ref, qn), u_trans_n(*ref, qn), tn, ta);
if(w == 0) continue;
kv_pushp(u_trans_t, *ta, &p);