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

Merge pull request #400 from chhylp123/hifiasm_dev_debug

Browse Source 
Hifiasm dev debug
This commit is contained in:
chhylp123
2023-02-16 00:09:39 -05:00
committed by GitHub
parent 1dd82e7729 286b897881
commit 026cf97d61
25 changed files with 4686 additions and 168 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+7
View File
@@ -1,3 +1,4 @@
#define __STDC_LIMIT_MACROS
#include <zlib.h>
#include <stdlib.h>
#include <stdio.h>
@@ -53,6 +54,8 @@ static ko_longopt_t long_options[] = {
{ "low-het", ko_no_argument, 339},
{ "s-base", ko_required_argument, 340},
{ "bin-only", ko_no_argument, 341},
{ "ul-round", ko_required_argument, 342},
{ "prt-raw", ko_no_argument, 343},
{ 0, 0, 0 }
};
@@ -262,6 +265,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->is_topo_trans = 1;
asm_opt->is_bub_trans = 1;
asm_opt->bin_only = 0;
asm_opt->ul_clean_round = 1;
asm_opt->prt_dbg_gfa = 0;
}
void destory_enzyme(enzyme* f)
@@ -794,6 +799,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
if(asm_opt->trans_base_rate_sec < 0) asm_opt->is_base_trans = 0;
}
else if (c == 341) asm_opt->bin_only = 1;
else if (c == 342) asm_opt->ul_clean_round = atol(opt.arg);
else if (c == 343) asm_opt->prt_dbg_gfa = 1;
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
+4 -1
View File
@@ -1,10 +1,11 @@
#ifndef __COMMAND_LINE_PARSER__
#define __COMMAND_LINE_PARSER__
#define __STDC_LIMIT_MACROS
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.18.5-r500"
#define HA_VERSION "0.18.6-r513"
#define VERBOSE 0
@@ -137,6 +138,8 @@ typedef struct {
uint8_t is_topo_trans;
uint8_t is_bub_trans;
uint8_t bin_only;
int32_t ul_clean_round;
int32_t prt_dbg_gfa;
} hifiasm_opt_t;
extern hifiasm_opt_t asm_opt;
Correct.cpp
+1461 -13
View File
File diff suppressed because it is too large Load Diff
Correct.h
+2
View File
@@ -1,5 +1,7 @@
#ifndef __CORRECT__
#define __CORRECT__
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include "Hash_Table.h"
#include "Levenshtein_distance.h"
Hash_Table.cpp
+1
View File
@@ -1,3 +1,4 @@
#define __STDC_LIMIT_MACROS
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
Levenshtein_distance.h
+2
View File
@@ -1,5 +1,7 @@
#ifndef __LEVENSHTEIN__
#define __LEVENSHTEIN__
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include "emmintrin.h"
#include "nmmintrin.h"
Output.h
+1
View File
@@ -1,6 +1,7 @@
#ifndef __OUTPUT__
#define __OUTPUT__
#define __STDC_LIMIT_MACROS
#include<stdint.h>
#include <string.h>
#include <stdlib.h>
Overlaps.cpp
+182 -56
View File
@@ -64,6 +64,7 @@ KSORT_INIT_GENERIC(uint32_t)
void reduce_hamming_error_adv(ma_ug_t *iug, asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub);
void print_vw_edge(asg_t *sg, uint32_t vid, uint32_t wid, const char *cmd);
typedef struct {
uint32_t d, tot, ma, p;
@@ -10101,13 +10102,29 @@ uint64_t *n_utg)
return C_bases/R_bases;
}
uint32_t cal_circle_ov(const ma_ug_t *ug, uint32_t uid, uint32_t rev, asg_t *sg)
{
uint32_t v, w, vx, wx, k;
v = w = (uid<<1)+(!!rev);
vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32));
wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32));
v = vx; w = wx;
if(sg) {
asg_arc_t *av; uint32_t nv;
av = asg_arc_a(sg, vx); nv = asg_arc_n(sg, vx);
for (k = 0; k < nv; k++) {
if(av[k].v == wx) return av[k].ol;
}
}
return 0;
}
void ma_ug_print2(const ma_ug_t *ug, All_reads *RNF, asg_t* read_g, const ma_sub_t *coverage_cut,
ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FILE *fp)
{
uint8_t* primary_flag = read_g?(uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t)):NULL;
uint32_t i, j, l, pc = read_g && coverage_cut && sources && ruIndex?1:0;
uint32_t i, j, l, pc = read_g && coverage_cut && sources && ruIndex?1:0, co;
char name[32];
for (i = 0; i < ug->u.n; ++i) { // the Segment lines in GFA
ma_utg_t *p = &ug->u.a[i];
@@ -10154,37 +10171,39 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FIL
uint32_t nu, u, v;
for (i = 0; i < ug->u.n; ++i) {
if(ug->u.a[i].m == 0) continue;
if(ug->u.a[i].circ)
{
fprintf(fp, "L\t%s%.6dc\t+\t%s%.6dc\t+\t%dM\tL1:i:%d\n",
prefix, i+1, prefix, i+1, 0, ug->u.a[i].len);
fprintf(fp, "L\t%s%.6dc\t-\t%s%.6dc\t-\t%dM\tL1:i:%d\n",
prefix, i+1, prefix, i+1, 0, ug->u.a[i].len);
if(ug->u.a[i].circ) {
co = cal_circle_ov(ug, i, 0, read_g);
fprintf(fp, "L\t%s%.6dc\t+\t%s%.6dc\t+\t%dM\tL1:i:%d\tL2:i:%u\n",
prefix, i+1, prefix, i+1, co, ((ug->u.a[i].len>=co)?(ug->u.a[i].len-co):0), 0);
co = cal_circle_ov(ug, i, 1, read_g);
fprintf(fp, "L\t%s%.6dc\t-\t%s%.6dc\t-\t%dM\tL1:i:%d\tL2:i:%u\n",
prefix, i+1, prefix, i+1, co, ((ug->u.a[i].len>=co)?(ug->u.a[i].len-co):0), 0);
} else {
u = i<<1;
au = asg_arc_a(ug->g, u);
nu = asg_arc_n(ug->g, u);
for (j = 0; j < nu; j++)
{
if(au[j].del) continue;
v = au[j].v;
fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\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]), 0/**au[j].ou**/);
}
u = (i<<1) + 1;
au = asg_arc_a(ug->g, u);
nu = asg_arc_n(ug->g, u);
for (j = 0; j < nu; j++)
{
if(au[j].del) continue;
v = au[j].v;
fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\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]), 0/**au[j].ou**/);
}
}
u = i<<1;
au = asg_arc_a(ug->g, u);
nu = asg_arc_n(ug->g, u);
for (j = 0; j < nu; j++)
{
if(au[j].del) continue;
v = au[j].v;
fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\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]), 0/**au[j].ou**/);
}
u = (i<<1) + 1;
au = asg_arc_a(ug->g, u);
nu = asg_arc_n(ug->g, u);
for (j = 0; j < nu; j++)
{
if(au[j].del) continue;
v = au[j].v;
fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\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]), 0/**au[j].ou**/);
}
}
}
free(primary_flag);
@@ -13494,12 +13513,16 @@ void hic_clean(asg_t* read_g)
void hic_clean_adv(asg_t *sg, ug_opt_t *uopt)
{
uint32_t i, k, m, z, v, w; ma_utg_t *u = NULL; uint32_t *ba, bn, n_vtx, beg, end, n0, n1;
uint32_t i, k, m, z, v, w, mk; ma_utg_t *u = NULL; uint32_t *ba, bn, n_vtx, beg, end, n0, n1; ma_utg_t *mz = NULL;
ma_ug_t *ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); n_vtx = ug->g->n_seq<<1; double bub_rate = 0.1;
uint8_t *bf = NULL; bubble_type *bub = gen_bubble_chain(sg, ug, uopt, &bf);
uint64_t tLen, vocc, socc, pocc; buf_t b; memset(&b, 0, sizeof(buf_t)); CALLOC(b.a, n_vtx);
REALLOC(bf, n_vtx); memset(bf, 0, sizeof((*bf))*n_vtx);
kvec_t(uint64_t) buf; kv_init(buf); n0 = n1 = 0;
for (i = 0; i < ug->g->n_seq; ++i) {
if(ug->g->seq[i].del) continue;
ug->g->seq[i].c = PRIMARY_LABLE;
}
for (i = 0; i < bub->b_ug->u.n; i++) {
u = &(bub->b_ug->u.a[i]);
@@ -13568,6 +13591,20 @@ void hic_clean_adv(asg_t *sg, ug_opt_t *uopt)
if((pocc+socc) >= (vocc*bub_rate)) continue;
pocc += socc;
asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0, NULL);
for (z = 0; z < b.b.n; z++) {
if(b.b.a[z]==v || b.b.a[z]==b.S.a[0]) continue;
// socc += ug->u.a[b.b.a[i]>>1].n;
if(ug->g->seq[b.b.a[z]>>1].del) continue;
if(ug->g->seq[b.b.a[z]>>1].c != ALTER_LABLE) continue;
mz = &(ug->u.a[b.b.a[z]>>1]);
if(mz->m == 0) continue;
for (mk = 0; mk < mz->n; mk++) asg_seq_del(sg, mz->a[mk]>>33);
asg_seq_del(ug->g, b.b.a[z]>>1);
if(ug->u.a[b.b.a[z]>>1].m) {
ug->u.a[b.b.a[z]>>1].m = ug->u.a[b.b.a[z]>>1].n = 0;
free(ug->u.a[b.b.a[z]>>1].a); ug->u.a[b.b.a[z]>>1].a = NULL;
}
}
// fprintf(stderr, "+utg%.6dl\tutg%.6dl\n", (int32_t)(v>>1)+1, (int32_t)(b.S.a[0]>>1)+1);
n0++;
}
@@ -13579,6 +13616,22 @@ void hic_clean_adv(asg_t *sg, ug_opt_t *uopt)
}
}
for (v = 0; v < ug->g->n_seq; ++v) {
if(ug->g->seq[v].del) continue;
if(ug->g->seq[v].c != ALTER_LABLE) continue;
mz = &(ug->u.a[v]);
if(mz->m == 0) continue;
for (k = 0; k < mz->n; k++) asg_seq_del(sg, mz->a[k]>>33);
asg_seq_del(ug->g, v);
if(ug->u.a[v].m) {
ug->u.a[v].m = ug->u.a[v].n = 0;
free(ug->u.a[v].a); ug->u.a[v].a = NULL;
}
}
asg_cleanup(ug->g);
tLen = get_bub_pop_max_dist_advance(ug->g, &b);
for (v = buf.n = 0; v < n_vtx; ++v) {
if(ug->g->seq[v>>1].del) continue;
@@ -13624,14 +13677,43 @@ void hic_clean_adv(asg_t *sg, ug_opt_t *uopt)
vocc += ug->u.a[b.b.a[k]>>1].n;
}
if((socc) >= (vocc*bub_rate)) continue;
asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0, NULL);
for (i = 0; i < b.b.n; i++) {
if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue;
// socc += ug->u.a[b.b.a[i]>>1].n;
if(ug->g->seq[b.b.a[i]>>1].del) continue;
if(ug->g->seq[b.b.a[i]>>1].c != ALTER_LABLE) continue;
mz = &(ug->u.a[b.b.a[i]>>1]);
if(mz->m == 0) continue;
for (mk = 0; mk < mz->n; mk++) asg_seq_del(sg, mz->a[mk]>>33);
asg_seq_del(ug->g, b.b.a[i]>>1);
if(ug->u.a[b.b.a[i]>>1].m) {
ug->u.a[b.b.a[i]>>1].m = ug->u.a[b.b.a[i]>>1].n = 0;
free(ug->u.a[b.b.a[i]>>1].a); ug->u.a[b.b.a[i]>>1].a = NULL;
}
}
// fprintf(stderr, "-utg%.6dl\tutg%.6dl\n", (int32_t)(v>>1)+1, (int32_t)(b.S.a[0]>>1)+1);
n1++;
}
}
}
filter_sg_by_ug(sg, ug, uopt);
// filter_sg_by_ug(sg, ug, uopt);
for (v = 0; v < ug->g->n_seq; ++v) {
if(ug->g->seq[v].del) continue;
if(ug->g->seq[v].c != ALTER_LABLE) continue;
mz = &(ug->u.a[v]);
if(mz->m == 0) continue;
for (k = 0; k < mz->n; k++) asg_seq_del(sg, mz->a[k]>>33);
asg_seq_del(ug->g, v);
if(ug->u.a[v].m) {
ug->u.a[v].m = ug->u.a[v].n = 0;
free(ug->u.a[v].a); ug->u.a[v].a = NULL;
}
}
asg_cleanup(ug->g);
asg_cleanup(sg);
free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
ma_ug_destroy(ug); free(bf); kv_destroy(buf);
destory_bubbles(bub); free(bub);
@@ -15627,7 +15709,6 @@ bub_label_t* b_mask_t)
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);
@@ -18754,10 +18835,11 @@ void filter_set_kug(uint8_t* trio_flag, asg_t *rg, uint8_t *rf, kvec_asg_arc_t_w
void output_trio_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, int is_bench, bub_label_t* b_mask_t)
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, int is_bench, long long gap_fuzz, ug_opt_t *opt, bub_label_t* b_mask_t)
{
reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, opt->ruIndex, NULL);
uint8_t *rf = NULL;
if(asm_opt.kpt_rate > 0) CALLOC(rf, sg->n_seq);
@@ -27602,7 +27684,7 @@ R_to_U* ruIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* new_rtg_edges,
void output_contig_graph_primary_pre(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, uint64_t bubble_dist, long long tipsLen,
R_to_U* ruIndex, int max_hang, int min_ovlp)
R_to_U* ruIndex, int max_hang, int min_ovlp, const ug_opt_t *uopt)
{
kvec_asg_arc_t_warp new_rtg_edges;
kv_init(new_rtg_edges.a);
@@ -27655,6 +27737,10 @@ R_to_U* ruIndex, int max_hang, int min_ovlp)
fclose(output_file);
}
///for debug
// graph_ovlp_binning(ug, sg, uopt);
// gen_hpc_re_t(ug);
free(gfa_name);
ma_ug_destroy(ug);
kv_destroy(new_rtg_edges.a);
@@ -31231,7 +31317,7 @@ int max_hang, int min_ovlp, uint32_t chainLenThres, long long gap_fuzz, bub_labe
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_backward(ug, sg, sources, coverage_cut, ruIndex, max_hang, min_ovlp, chainLenThres, beg_idx, occ, &bub, b_mask_t);
/**
if((!no_trio_recover) && (ha_opt_triobin(&asm_opt)))
{
ma_ug_destroy(ug); ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
@@ -31239,6 +31325,7 @@ int max_hang, int min_ovlp, uint32_t chainLenThres, long long gap_fuzz, bub_labe
// rescue_missing_hap_ovlp(ug, sg, sources, coverage_cut, max_hang, min_ovlp, &bub, gap_fuzz);
reduce_hamming_error_adv(ug, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, ruIndex, &bub);
}
**/
destory_bubbles(&bub);
destory_hap_cov_t(&cov);
@@ -33582,10 +33669,10 @@ ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t,
}
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)
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, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, uint64_t deep_clean)
{
ul_renew_t nopt; memset(&nopt, 0, sizeof(nopt));
nopt.src = sources; nopt.r_src = reverse_sources; nopt.ruIndex = ruIndex;
@@ -33593,7 +33680,7 @@ int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file)
nopt.cov = coverage_cut; nopt.b_mask_t = b_mask_t;
nopt.max_hang = max_hang_length; nopt.mini_ovlp = mini_overlap_length;
ul_realignment_gfa(uopt, *sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
asm_opt.max_short_tip, asm_opt.max_short_ul_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file, &nopt);;
asm_opt.max_short_tip, asm_opt.max_short_ul_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file, &nopt, bin_file, free_uld, is_bridg, deep_clean);
// 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);
@@ -33607,6 +33694,40 @@ int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file)
// post_rescue(uopt, *sg, (*sources), (*reverse_sources), ruIndex, b_mask_t, 0);
}
void gradually_renew_g(ma_hit_t_alloc **src, ma_hit_t_alloc **rev_src, long long *n_read,
uint64_t **readLen, ma_sub_t **cov, 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, int64_t gap_fuzz, int64_t min_dp,
bub_label_t *b_mask_t, uint32_t is_trio, int32_t ul_aln_round, char *o_file, const char *bin_file)
{
int32_t k, strl = strlen(bin_file)+1, kt, cl, sl; char *id = NULL;
renew_g(src, rev_src, n_read, readLen, cov, 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,
is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0, 1, ((is_trio)?(0):(1))/**1**/);
gen_ug_opt_t(uopt, *src, *rev_src, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, *readLen,
*cov, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, b_mask_t);
ug_ext_gfa(uopt, *sg, ug_ext_len);
/**if(!ha_opt_triobin(&asm_opt))**/ hic_clean_adv(*sg, uopt);
cl = strl+1; MALLOC(id, cl);
for (k = 1; k < ul_aln_round; k++) {
for(kt = k, sl = strl+1; kt > 0; kt/=10) sl++;
if(cl < sl) {
cl = sl; REALLOC(id, cl);
}
sprintf(id, "%s%d", bin_file, k);
renew_g(src, rev_src, n_read, readLen, cov, 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,
is_trio, o_file, id, ((k+1)==ul_aln_round)?1:0, 0, 0);
gen_ug_opt_t(uopt, *src, *rev_src, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, *readLen,
*cov, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, b_mask_t);
ug_ext_gfa(uopt, *sg, ug_ext_len);
/**if(!ha_opt_triobin(&asm_opt))**/ hic_clean_adv(*sg, uopt);
}
free(id);
}
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,
@@ -33711,6 +33832,20 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
// set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length);
}
if(asm_opt.ar) {
gradually_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, gap_fuzz, min_dp, &b_mask_t,
ha_opt_triobin(&asm_opt), asm_opt.ul_clean_round, o_file, "re");
} else {
ug_ext_gfa(&uopt, sg, ug_ext_len);
if(!ha_opt_triobin(&asm_opt)) {
// output_unitig_graph(sg, coverage_cut, "pre_clean", sources, ruIndex, max_hang_length, mini_overlap_length);
hic_clean_adv(sg, &uopt);
}
}
/**
if(asm_opt.ar) {
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,
@@ -33718,11 +33853,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
///make sure uopt has been updated; not necessary
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);
// 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);
/**
asg_cut_tip(sg, asm_opt.max_short_tip);
@@ -33853,15 +33985,9 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
// flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL;
flat_soma_v(sg, sources, ruIndex);
**/
ug_ext_gfa(&uopt, sg, ug_ext_len);
if(!ha_opt_triobin(&asm_opt)) {
// output_unitig_graph(sg, coverage_cut, "pre_clean", sources, ruIndex, max_hang_length, mini_overlap_length);
hic_clean_adv(sg, &uopt);
}
output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources,
asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length);
asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length, &uopt);
/**
if (asm_opt.flag & HA_F_VERBOSE_GFA)
{
@@ -33885,7 +34011,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
{
if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION;
output_trio_graph(sg, coverage_cut, o_file, sources, reverse_sources, (asm_opt.max_short_tip*2),
0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 0, &b_mask_t);
0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 0, gap_fuzz, &uopt, &b_mask_t);
}
else if(ha_opt_hic(&asm_opt))
{
Overlaps.h
+45 -3
View File
@@ -1,5 +1,7 @@
#ifndef __OVERLAPS__
#define __OVERLAPS__
#define __STDC_LIMIT_MACROS
#include <stdio.h>
#include <stdint.h>
#include "kvec.h"
@@ -32,7 +34,6 @@
// #define PRIMARY_LABLE 1
// #define ALTER_LABLE 2
// #define HAP_LABLE 4
#define HA_RE_UL_ID "re"
#define ug_ext_len 75000
#define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32))
@@ -68,6 +69,22 @@ typedef struct {
size_t n, m;
} kv_ul_ov_t;
typedef struct {
uint32_t tn, rn, el;
uint32_t qs, qe, ts, te;
uint8_t dir:5, pe:1, full:1, rev:1;
} emask_t;
typedef struct {
emask_t *a;
size_t n, m;
} kv_emask_t;
typedef struct {
kv_emask_t *a;
uint32_t n;
} idx_emask_t;
typedef struct {
///off: start idx in mg128_t * a[];
///cnt: how many eles in this chain
@@ -119,8 +136,6 @@ void ma_hit_sort_qns(ma_hit_t *a, long long n);
int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources,
char* output_file_name);
typedef struct {
uint32_t s:31, del:1, e;
uint8_t c;
@@ -240,6 +255,7 @@ void print_gfa(asg_t *g);
typedef struct { size_t n, m; uint64_t *a; } asg64_v;
typedef struct { size_t n, m; uint32_t *a; } asg32_v;
typedef struct { size_t n, m; ma_utg_t *a;} ma_utg_v;
typedef struct { asg64_v idx; kv_ul_ov_t srt;} mask_ul_ov_t;
typedef struct {
ma_utg_v u;
@@ -295,9 +311,33 @@ typedef struct {
hmap_t *mm;
} hpc_t;
typedef struct {
uint32_t s, e;
uint8_t k;
} hpc_ss_t;
typedef struct {
uint32_t s, e;
} hpc_idx_t;
typedef struct {
size_t n, m;
hpc_ss_t *a;
hpc_idx_t *idx;
uint64_t idx_n;
} hpc_re_t;
#define hpc_len(x, id) ((x).hg->u.a[(id)].len>>1)
#define hpc_str(x, id, rev) (((x).hg->u.a[(id)].s)+((rev)?((x).hg->u.a[(id)].len>>1):(0)))
typedef struct {
uint32_t n;
uint8_t *a;
} bit_mask_t;
#define set_bit_mask_t(x, i) ((x).a[(i)>>3]|=(((uint8_t)1)<<((i)&((uint32_t)7))))
#define get_bit_mask_t(x, i) ((x).a[(i)>>3]&(((uint8_t)1)<<((i)&((uint32_t)7))))
typedef struct {
ma_ug_t *ug;
hpc_t *hpc_g;
@@ -1163,6 +1203,8 @@ void set_reverse_overlap(ma_hit_t* dest, ma_hit_t* source);
// int* b_low_cov, int* b_high_cov, double m_rate);
void ma_hit_contained_advance(ma_hit_t_alloc* sources, long long n_read, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp);
void hic_clean_adv(asg_t *sg, ug_opt_t *uopt);
void update_ug_ou(ma_ug_t *ug, asg_t *sg);
#define JUNK_COV 5
#define DISCARD_RATE 0.8
POA.h
+2
View File
@@ -1,5 +1,7 @@
#ifndef __POA_PARSER__
#define __POA_PARSER__
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include "Hash_Table.h"
#include "Process_Read.h"
Process_Read.cpp
+5
View File
@@ -764,6 +764,11 @@ void destory_all_ul_t(all_ul_t *x) {
for (i = 0; i < x->nid.n; i++) free(x->nid.a[i].a);
free(x->nid.a);
free(x->ridx.idx.a); free(x->ridx.occ.a);
// if(x->mm) {
// for (i = 0; i < x->mm->n; i++) free(x->mm->a[i].a);
// free(x->mm->a); free(x->mm); x->mm = NULL;
// }
}
void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn, uint64_t nn_offset)
Process_Read.h
+2
View File
@@ -1,6 +1,7 @@
#ifndef __READ__
#define __READ__
#define __STDC_LIMIT_MACROS
#include<stdint.h>
#include <string.h>
#include <stdlib.h>
@@ -198,6 +199,7 @@ typedef struct
ul_vec_t *a;
size_t n, m;
All_reads *hR;
// idx_emask_t *mm;
// uint32_t mm;
} all_ul_t;
Purge_Dups.h
+2
View File
@@ -1,5 +1,7 @@
#ifndef __PURGEDUPS__
#define __PURGEDUPS__
#define __STDC_LIMIT_MACROS
#include <stdio.h>
#include <stdint.h>
#include "kvec.h"
anchor.cpp
+20
View File
@@ -1806,3 +1806,23 @@ int64_t ug_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint6
return 0;
}
}
int64_t ug_map_lchain_simple(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres,
int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate,
uint32_t gen_off, double mcopy_rate, uint32_t mcopy_khit_cut, uint32_t is_hpc, ha_mzl_t *res, uint64_t res_n, ha_mzl_t *idx, uint64_t idx_n, uint64_t mzl_cutoff, uint64_t chain_cutoff)
{
int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip;
set_lchain_dp_op(is_accurate, mz_k, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check);
if((!overlap_list) || (!cl)) {
ab->mz.n = 0;
mz2_ha_sketch(rs, rl, mz_w, mz_k, rid, is_hpc, &ab->mz, NULL, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL);
if(res) memcpy(res, ab->mz.a, ab->mz.n * (sizeof((*(ab->mz.a)))));
return ab->mz.n;
} else {
sp->n = 0;
minimizers_qgen_input(ab, rid, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, NULL, uref, dbg_ct, sp, high_occ, low_occ, res, res_n, idx, idx_n, mzl_cutoff, chain_cutoff, NULL);
// lchain_qgen_mcopy_input(cl, overlap_list, rid, rl, NULL, uref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, bw_thres_sec, quick_check, gen_off, mcopy_rate, mcopy_khit_cut, sp);
lchain_qgen_mcopy(cl, overlap_list, rid, rl, NULL, uref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off, mcopy_rate, chain_cutoff, mcopy_khit_cut, sp);
return 0;
}
}
gfa_ut.cpp
+813 -69
View File
File diff suppressed because it is too large Load Diff
gfa_ut.h
+1 -1
View File
@@ -27,7 +27,7 @@ void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float o
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,
double max_ovlp_drop_ratio, int64_t max_tip, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, ul_renew_t *ropt);
double max_ovlp_drop_ratio, int64_t max_tip, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, ul_renew_t *ropt, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, uint64_t deep_clean);
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);
hic.h
+2
View File
@@ -1,5 +1,7 @@
#ifndef __HIC__
#define __HIC__
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include "Overlaps.h"
horder.cpp
+7
View File
@@ -838,6 +838,13 @@ void horder_clean_sg_by_utg(asg_t *sg, ma_ug_t *ug)
if(av[k].del) continue;
w = av[k].v;
vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32));
wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32));
asg_arc_unique_del(sg, vx, wx, 0); asg_arc_unique_del(sg, wx^1, vx^1, 0);
}
if(u->circ) {
v = w = i<<1;
vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32));
wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32));
asg_arc_unique_del(sg, vx, wx, 0); asg_arc_unique_del(sg, wx^1, vx^1, 0);
horder.h
+2
View File
@@ -1,5 +1,7 @@
#ifndef __HORDER__
#define __HORDER__
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include "hic.h"
inter.cpp
+2112 -24
View File
File diff suppressed because it is too large Load Diff
inter.h
+7 -1
View File
@@ -108,7 +108,7 @@ void ul_resolve(ma_ug_t *ug, const asg_t *rg, const ug_opt_t *uopt, int hap_n);
void ul_load(const ug_opt_t *uopt);
uint64_t* get_hifi2ul_list(all_ul_t *x, uint64_t hid, uint64_t* a_n);
uint64_t ul_refine_alignment(const ug_opt_t *uopt, asg_t *sg);
ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_cache);
ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_cache, const char *bin_file);
int32_t write_all_ul_t(all_ul_t *x, char* file_name, ma_ug_t *ug);
int32_t load_all_ul_t(all_ul_t *x, char* file_name, All_reads *hR, ma_ug_t *ug);
uint32_t ugl_cover_check(uint64_t is, uint64_t ie, ma_utg_t *u);
@@ -118,5 +118,11 @@ void update_ug_arch_ul_mul(ma_ug_t *ug);
void print_ul_alignment(ma_ug_t *ug, all_ul_t *aln, uint32_t id, const char* cmd);
void clear_all_ul_t(all_ul_t *x);
void trans_base_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res, bubble_type *bub);
hpc_re_t *gen_hpc_re_t(ma_ug_t *ug);
idx_emask_t* graph_ovlp_binning(ma_ug_t *ug, asg_t *sg, const ug_opt_t *uopt);
uint32_t gen_src_shared_interval_simple(uint32_t src, ma_ug_t *ug, kv_ul_ov_t *res);
uint64_t check_ul_ov_t_consist(ul_ov_t *x, ul_ov_t *y, int64_t ql, int64_t tl, double diff);
uint32_t infer_se(uint32_t qs, uint32_t qe, uint32_t ts, uint32_t te, uint32_t rev,
uint32_t rqs, uint32_t rqe, uint32_t *rts, uint32_t *rte);
#endif
ksw2.h
+1
View File
@@ -1,6 +1,7 @@
#ifndef KSW2_H_
#define KSW2_H_
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#define KSW_NEG_INF -0x40000000
kthread.cpp
+1
View File
@@ -1,3 +1,4 @@
#define __STDC_LIMIT_MACROS
#include <pthread.h>
#include <stdlib.h>
#include <limits.h>
rcut.h
+2
View File
@@ -1,5 +1,7 @@
#ifndef __RCUT__
#define __RCUT__
#define __STDC_LIMIT_MACROS
#include <stdio.h>
#include <stdint.h>
#include "kvec.h"
tovlp.h
+2
View File
@@ -1,5 +1,7 @@
#ifndef __TOVLP__
#define __TOVLP__
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include "Overlaps.h"