unsucessful hpc mark

This commit is contained in:
chhylp123
2026-05-15 19:02:16 -04:00
parent 7884b5ad88
commit bca2394b69
12 changed files with 2679 additions and 216 deletions
+20 -12
View File
@@ -1023,10 +1023,10 @@ void ha_ec(int64_t round, int num_pround, int des_idx, uint64_t *tot_b, uint64_t
if(round == asm_opt.number_of_round-1 && asm_opt.is_dbg_het_cnt) CALLOC(het_cnt, R_INF.total_reads);
if (r_out) {
write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, 0);
if((asm_opt.flag & HA_F_VERBOSE_GFA) && (asm_opt.bin_only == 1)) exit(1);///just for debug
}
if (w_tmp) tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, round, asm_opt.number_of_round, 0);
if (w_tmp) tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, round, asm_opt.number_of_round, 0, 0);
// Output_corrected_fastq();
@@ -1112,7 +1112,7 @@ void ha_overlap_and_correct(int round)
// fprintf(stderr, "[M::%s::%.3f] ==> chaining\n", __func__, yak_realtime_0()-tt0);
// exit(1);
if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, 0);
ha_pt_destroy(ha_idx);
ha_idx = NULL;
@@ -1727,6 +1727,7 @@ void Output_PAF()
fclose(output_file);
fprintf(stderr, "PAF has been written.\n");
exit(1);
}
@@ -1952,12 +1953,12 @@ void ha_overlap_final(void)
asm_opt.het_cov = het_cov;
}
void ha_ec_ff(int renew_idx)
void ha_ec_ff(int renew_idx, int8_t pre_load_idx, uint64_t w_tmp)
{
int hom_cov, het_cov;
ha_flt_tab_hp = ha_idx_hp = NULL;
if(ha_idx && renew_idx) {
if((ha_idx) && (renew_idx) && (!pre_load_idx)) {
ha_pt_destroy(ha_idx); ha_idx = NULL;
}
@@ -1966,6 +1967,8 @@ void ha_ec_ff(int renew_idx)
asm_opt.hom_cov = hom_cov; asm_opt.het_cov = het_cov;
}
if (w_tmp) tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, asm_opt.number_of_round, asm_opt.number_of_round, 0, 1);
cal_ov_r(asm_opt.thread_num, R_INF.total_reads, renew_idx);
if(asm_opt.write_pos_idx) {
@@ -2076,7 +2079,7 @@ int ha_assemble(void)
// 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, r0 = -1, hom_cov = -1, ovlp_loaded = 0; uint64_t tot_b, tot_e;
int r, r0 = -1, hom_cov = -1, ovlp_loaded = 0; uint64_t tot_b, tot_e; int8_t pre_load_ff = 0;
if ((asm_opt.load_index_from_disk) && (asm_opt.dbg_ec_rr < 0) && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
ovlp_loaded = 1;
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
@@ -2093,13 +2096,13 @@ int ha_assemble(void)
}
if (!ovlp_loaded) {
ha_flt_tab = ha_idx = NULL;
if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name), load_ct_index(&ha_ct_table, asm_opt.output_file_name);
if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, 0), load_ct_index(&ha_ct_table, asm_opt.output_file_name);
r = ha_idx?asm_opt.number_of_round-1:0;
if((!ha_idx) && (asm_opt.restart)) {
r = asm_opt.number_of_round - 1;
r = asm_opt.number_of_round/** - 1**/;
if(asm_opt.dbg_ec_rr >= 0) r = asm_opt.dbg_ec_rr;
for (; r >= 0; --r) {
if(tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, r, asm_opt.number_of_round, 1)) {
if(tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, r, asm_opt.number_of_round, 1, ((r==asm_opt.number_of_round)?(1):(0)))) {
load_ct_index(&ha_ct_table, asm_opt.output_file_name); r0 = r;
break;
}
@@ -2108,7 +2111,11 @@ int ha_assemble(void)
fprintf(stderr, "[E::%s] no matching debug error-correction bins found\n", __func__);
exit(1);
}
if(r < 0) r = 0;
if(r < 0) {
r = 0;
} else if(r == asm_opt.number_of_round) {
pre_load_ff = 1;
}
}
// construct hash table for high occurrence k-mers
@@ -2138,8 +2145,9 @@ int ha_assemble(void)
// overlap between corrected reads
ha_opt_reset_to_round(&asm_opt, asm_opt.number_of_round);
// ha_overlap_final();
ha_ec_ff(1/**0**/);
ha_ec_ff(1/**0**/, pre_load_ff, ((r > r0) && (asm_opt.restart))?1:0);
fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb());
if(asm_opt.dbg_ec_rr >= 0) exit(1);
// fprintf(stderr, "\n[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb());
// ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads);
if(!(asm_opt.write_pos_idx)) {
@@ -2193,7 +2201,7 @@ int ha_assemble_pair(void)
// Output_corrected_reads(); exit(0);
ha_flt_tab = ha_idx = NULL; r = asm_opt.number_of_round - 1;
if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name), load_ct_index(&ha_ct_table, asm_opt.output_file_name);
if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, 0), load_ct_index(&ha_ct_table, asm_opt.output_file_name);
// construct hash table for high occurrence k-mers
if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL) {
+127 -3
View File
@@ -12,6 +12,11 @@
KSEQ_INIT(gzFile, gzread)
#define DEFAULT_OUTPUT "hifiasm.asm"
#define RER_H_HiFi 0.06
#define RER_H_ONT 0.08
#define RER_N_HiFi 0.03
#define RER_N_ONT 0.05
#define MIN_R_COV 3
hifiasm_opt_t asm_opt;
@@ -94,7 +99,14 @@ static ko_longopt_t long_options[] = {
{ "hyb-syn", ko_required_argument, 376},
{ "simd-m", ko_required_argument, 377},
{ "del-hf", ko_no_argument, 378},
{ "dbg-rr", ko_required_argument, 379},
{ "dbg-rr", ko_required_argument, 379},///start from 0/1/2 step of correction
{ "re-aln", ko_required_argument, 380},
{ "syn", ko_required_argument, 381},
{ "rec", ko_required_argument, 382},
{ "ref", ko_required_argument, 383},
{ "h_rec", ko_required_argument, 384},
{ "h_ref", ko_required_argument, 385},
{ "ret", ko_no_argument, 386},
// { "path-round", ko_required_argument, 348},
{ 0, 0, 0 }
};
@@ -106,6 +118,7 @@ double Get_T(void)
return t.tv_sec+t.tv_usec/1000000.0;
}
void Print_H(hifiasm_opt_t* asm_opt)
{
fprintf(stderr, "Usage: hifiasm [options] <in_1.fq> <in_2.fq> <...>\n");
@@ -140,7 +153,34 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " discard overlaps supported by <INT minimizers [%ld]\n", asm_opt->chn_occ);
fprintf(stderr, " --ec-only error correction only; disable overlapping and assembly\n");
fprintf(stderr, " --simd-m use SIMD acceleration when supported: AVX-512 (2), AVX2 (1), or non-SIMD (0)\n");
fprintf(stderr, " --re-aln keep raw reads for error correction [%d]\n", asm_opt->realn_raw);
fprintf(stderr, " --syn rescue overlaps for error correction [%d]\n", asm_opt->post_syn);
fprintf(stderr, " Recurrent sequencing-error filtering:\n");
fprintf(stderr, " --ret enable recurrent sequencing-error filtering; disabled by default\n");
fprintf(stderr, " thresholds are controlled by --rec/--ref and --h_rec/--h_ref\n");
fprintf(stderr, " --rec/--ref and --h_rec/--h_ref are ignored unless --ret is set\n");
fprintf(stderr, " --rec INT\n");
fprintf(stderr, " discard candidate variants supported by <INT reads [%ld]\n",
asm_opt->recurrent_err_normal_min);
fprintf(stderr, " --ref FLOAT\n");
fprintf(stderr, " discard candidate variants supported by <FLOAT*coverage reads; -1 to disable\n");
fprintf(stderr, " coverage = min(local read coverage, homozygous coverage);\n");
fprintf(stderr, " default: %.3g for ONT, %.3g for HiFi\n",
RER_N_ONT, RER_N_HiFi);
fprintf(stderr, " --h_rec INT\n");
fprintf(stderr, " discard candidate variants near homopolymers if supported by <INT reads [%ld]\n",
asm_opt->recurrent_err_hpc_min);
fprintf(stderr, " --h_ref FLOAT\n");
fprintf(stderr, " discard candidate variants near homopolymers if supported by <FLOAT*coverage reads; -1 to disable\n");
fprintf(stderr, " coverage = min(local read coverage, homozygous coverage);\n");
fprintf(stderr, " default: %.3g for ONT, %.3g for HiFi\n",
RER_H_ONT, RER_H_HiFi);
fprintf(stderr, " Assembly:\n");
fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round);
fprintf(stderr, " -m INT pop bubbles of <INT in size in contig graphs [%lld]\n", asm_opt->large_pop_bubble_size);
@@ -365,8 +405,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->dp_e = 0.0025;
asm_opt->hg_size = -1;
asm_opt->kpt_rate = -1;
asm_opt->infor_cov = 3;
asm_opt->s_hap_cov = 3;
asm_opt->infor_cov = MIN_R_COV;
asm_opt->s_hap_cov = MIN_R_COV;
asm_opt->ul_error_rate = 0.2/**0.15**/;
asm_opt->ul_error_rate_low = 0.1;
asm_opt->ul_error_rate_hpc = 0.2;
@@ -447,6 +487,22 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->del_hf = 0;
asm_opt->dbg_ec_rr = -1;
asm_opt->realn_raw = 0;
asm_opt->post_syn = 1;
asm_opt->recurrent_err_normal_min = MIN_R_COV;
asm_opt->recurrent_err_normal_min_set = -1;
asm_opt->recurrent_err_normal_rat = RER_N_HiFi;
asm_opt->recurrent_err_normal_rat_set = -1;
asm_opt->recurrent_err_hpc_min = MIN_R_COV;
asm_opt->recurrent_err_hpc_min_set = -1;
asm_opt->recurrent_err_hpc_rat = RER_H_HiFi;
asm_opt->recurrent_err_hpc_rat_set = -1;
asm_opt->recurrent_err_test = 0;
}
void destory_enzyme(enzyme* f)
@@ -837,6 +893,26 @@ int check_option(hifiasm_opt_t* asm_opt)
return 0;
}
if((asm_opt->realn_raw != 0) && (asm_opt->realn_raw != 1)) {
fprintf(stderr, "[ERROR] [--re-aln] must be 0/1\n");
return 0;
}
if((asm_opt->post_syn != 0) && (asm_opt->post_syn != 1)) {
fprintf(stderr, "[ERROR] [--syn] must be 0/1\n");
return 0;
}
if(asm_opt->recurrent_err_normal_rat > 1.0) {
fprintf(stderr, "[ERROR] [--ref] must be >= 0 && <= 1.0; -1 to disable\n");
return 0;
}
if(asm_opt->recurrent_err_hpc_rat > 1.0) {
fprintf(stderr, "[ERROR] [--h_ref] must be >= 0 && <= 1.0; -1 to disable\n");
return 0;
}
return 1;
}
@@ -1126,6 +1202,24 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
asm_opt->del_hf = 1;
} else if (c == 379) {
asm_opt->dbg_ec_rr = atoi(opt.arg);
} else if (c == 380) {
asm_opt->realn_raw = atoi(opt.arg);
} else if (c == 381) {
asm_opt->post_syn = atoi(opt.arg);
} else if (c == 382) {
asm_opt->recurrent_err_normal_min = atoi(opt.arg);
asm_opt->recurrent_err_normal_min_set = 1;
} else if (c == 383) {
asm_opt->recurrent_err_normal_rat = atof(opt.arg);
asm_opt->recurrent_err_normal_rat_set = 1;
} else if (c == 384) {
asm_opt->recurrent_err_hpc_min = atoi(opt.arg);
asm_opt->recurrent_err_hpc_min_set = 1;
} else if (c == 385) {
asm_opt->recurrent_err_hpc_rat = atof(opt.arg);
asm_opt->recurrent_err_hpc_rat_set = 1;
} else if (c == 386) {
asm_opt->recurrent_err_test = 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);
}
@@ -1167,5 +1261,35 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
asm_opt->rl_cut = -1; asm_opt->sc_cut = 1;
}
if(asm_opt->recurrent_err_normal_min_set != 1) {
asm_opt->recurrent_err_normal_min = MIN_R_COV;
}
if(asm_opt->recurrent_err_normal_rat_set != 1) {
asm_opt->recurrent_err_normal_rat = ((asm_opt->is_ont)?(RER_N_ONT):(RER_N_HiFi));
}
if(asm_opt->recurrent_err_hpc_min_set != 1) {
asm_opt->recurrent_err_hpc_min = MIN_R_COV;
}
if(asm_opt->recurrent_err_hpc_rat_set != 1) {
asm_opt->recurrent_err_hpc_rat = ((asm_opt->is_ont)?(RER_H_ONT):(RER_H_HiFi));
}
if(asm_opt->recurrent_err_normal_min < 0) asm_opt->recurrent_err_normal_min = -1;
if(asm_opt->recurrent_err_normal_rat < 0) asm_opt->recurrent_err_normal_rat = -1;
if(asm_opt->recurrent_err_hpc_min < 0) asm_opt->recurrent_err_hpc_min = -1;
if(asm_opt->recurrent_err_hpc_rat < 0) asm_opt->recurrent_err_hpc_rat = -1;
if(asm_opt->recurrent_err_test == 0) {
asm_opt->recurrent_err_normal_min = -1;
asm_opt->recurrent_err_normal_min_set = -1;
asm_opt->recurrent_err_normal_rat = -1;
asm_opt->recurrent_err_normal_rat_set = -1;
asm_opt->recurrent_err_hpc_min = -1;
asm_opt->recurrent_err_hpc_min_set = -1;
asm_opt->recurrent_err_hpc_rat = -1;
asm_opt->recurrent_err_hpc_rat_set = -1;
}
return check_option(asm_opt);
}
+19 -1
View File
@@ -5,7 +5,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.25.1-r920"
#define HA_VERSION "0.25.1-r933"
#define VERBOSE 0
@@ -206,6 +206,24 @@ typedef struct {
int8_t del_hf;
int64_t dbg_ec_rr;
int8_t realn_raw;
int8_t post_syn;
int64_t recurrent_err_normal_min;
int8_t recurrent_err_normal_min_set;
double recurrent_err_normal_rat;
int8_t recurrent_err_normal_rat_set;
int64_t recurrent_err_hpc_min;
int8_t recurrent_err_hpc_min_set;
double recurrent_err_hpc_rat;
int8_t recurrent_err_hpc_rat_set;
int8_t recurrent_err_test;
} hifiasm_opt_t;
extern hifiasm_opt_t asm_opt;
+1506 -108
View File
File diff suppressed because it is too large Load Diff
+1 -1
View File
@@ -1458,7 +1458,7 @@ void gen_hc_r_alin_nec_adv(gen_hc_aln_t *ez);
uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
uint64_t gen_hc_r_alin_self(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, asg16_v* scc);
void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, asg16_v *qc0, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8, uint8_t lindel, uint64_t tcut, uint64_t site_sc, int64_t h0_w, asg32_v *b32,
int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate, overlap_region *rchn, int64_t rref_len);
int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate, overlap_region *rchn, int64_t rref_len, int64_t flag_hf_ov_cut, int64_t min_re_cut, double min_re_rt, int64_t min_hpc_re_cut, double min_hpc_re_rt, int8_t is_hpc_flt);
void set_exact_exz(bit_extz_t *exz, int64_t qs, int64_t qe, int64_t ts, int64_t te);
void push_alnw(overlap_region *aux_o, bit_extz_t *exz);
void cal_exz_global(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez);
+148 -17
View File
@@ -178,6 +178,12 @@ typedef struct {
int64_t min_sc, penalty, max_drop;
} telo_end_pip_t;
typedef struct {
ma_hit_t_alloc* src;
ma_hit_t_alloc* rec;
int64_t n_thread, n_a;
} rclean_aux;
///this value has been updated at the first line of build_string_graph_without_clean
long long min_thres;
@@ -188,6 +194,24 @@ kv_u_trans_t *get_utg_ovlp(ma_ug_t **ug, asg_t* read_g, ma_hit_t_alloc* sources,
R_to_U* ruIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t, uint8_t* r_het);
void delete_useless_nodes(ma_ug_t **ug);
static void reset_ma_hit_t_0(void *data, long i, int tid) // callback for kt_for()
{
rclean_aux *p = (rclean_aux *)data; uint64_t k;
for (k = 0; k < p->src[i].length; k++) {
p->src[i].buffer[k].del = 0;
}
for (k = 0; k < p->rec[i].length; k++) {
p->rec[i].buffer[k].del = 0;
}
}
void inital_ovlap(ma_hit_t_alloc *src, ma_hit_t_alloc *rec, int64_t n_a, int64_t n_thre)
{
rclean_aux p; memset(&p, 0, sizeof(p));
p.n_thread = n_thre; p.src = src; p.rec = rec; p.n_a = n_a;
kt_for(p.n_thread, reset_ma_hit_t_0, &p, p.n_a);
}
static void mark_telo_ends(void *data, long i, int tid) // callback for kt_for()
{
telo_end_pip_t *sl = (telo_end_pip_t *)data;
@@ -23758,7 +23782,7 @@ void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sou
}
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, all_ul_t *ul_r_inf);
char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_ul_t *ul_r_inf, uint8_t **cmk);
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);
@@ -23776,7 +23800,7 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s
}
}
if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_debug_graph(NULL, NULL, NULL, output_file_name, NULL, NULL, NULL)) {
if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_debug_graph(NULL, NULL, NULL, output_file_name, NULL, NULL, NULL, NULL)) {
(*sources) = NULL;
(*reverse_sources) = NULL;
free(gfa_name);
@@ -32191,6 +32215,42 @@ int load_coverage_cut(ma_sub_t** coverage_cut, char* read_file_name)
}
int write_cmk(uint8_t* cmk, char* read_file_name, uint64_t n_read)
{
char* index_name = (char*)malloc(strlen(read_file_name)+15);
sprintf(index_name, "%s.bin", read_file_name);
FILE* fp = fopen(index_name, "w");
fwrite(&n_read, sizeof(n_read), 1, fp);
fwrite(cmk, sizeof((*(cmk))), n_read, fp);
free(index_name);
fflush(fp);
fclose(fp);
return 1;
}
int load_cmk(uint8_t** cmk, char* read_file_name)
{
char* index_name = (char*)malloc(strlen(read_file_name)+15);
sprintf(index_name, "%s.bin", read_file_name);
FILE* fp = fopen(index_name, "r");
if(!fp)
{
return 0;
}
int f_flag = 0;
uint64_t n_read;
f_flag += fread(&n_read, sizeof(n_read), 1, fp);
(*cmk) = (uint8_t*)malloc(sizeof((*(*(cmk))))*n_read);
fread((*cmk), sizeof((*((*cmk)))), n_read, fp);
free(index_name);
fflush(fp);
fclose(fp);
return 1;
}
int write_coverage_cut(ma_sub_t* coverage_cut, char* read_file_name, uint64_t n_read)
{
char* index_name = (char*)malloc(strlen(read_file_name)+15);
@@ -32441,7 +32501,7 @@ int load_asg_t(asg_t **sg, char* read_file_name)
}
int write_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, all_ul_t *ul_r_inf)
char* output_file_name, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, all_ul_t *ul_r_inf, uint8_t *cmk)
{
char* gfa_name = (char*)malloc(strlen(output_file_name)+55);
@@ -32465,6 +32525,10 @@ char* output_file_name, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, all_ul
sprintf(gfa_name, "%s.all.debug.ul.rinfor", output_file_name);
write_all_ul_t(ul_r_inf, gfa_name, NULL);
}
if(cmk) {
sprintf(gfa_name, "%s.all.debug.cmk", output_file_name);
write_cmk(cmk, gfa_name, R_INF.total_reads);
}
free(gfa_name);
fprintf(stderr, "debug_graph has been written.\n");
return 1;
@@ -32472,7 +32536,7 @@ char* output_file_name, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, all_ul
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, all_ul_t *ul_r_inf)
char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_ul_t *ul_r_inf, uint8_t **cmk)
{
FILE* fp = NULL;
char* gfa_name = (char*)malloc(strlen(output_file_name)+55);
@@ -32494,6 +32558,11 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u
sprintf(gfa_name, "%s.all.debug.ul.rinfor.ul.ovlp.bin", output_file_name);
fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp);
}
if(cmk) {
sprintf(gfa_name, "%s.all.debug.cmk.bin", output_file_name);
fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp);
}
if((sources == NULL) || (reverse_sources == NULL) || (ruIndex == NULL))
{
@@ -32567,6 +32636,11 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u
if(!load_all_ul_t(ul_r_inf, gfa_name, &R_INF, NULL)) return 0;
}
if(cmk) {
sprintf(gfa_name, "%s.all.debug.cmk", output_file_name);
if(!load_cmk(cmk, gfa_name)) return 0;
}
R_INF.paf = (*sources); R_INF.reverse_paf = (*reverse_sources);
return 1;
@@ -39767,19 +39841,24 @@ ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t,
// prt_dbg_rid_ovlp(src, -1, (char*)"c85c2e91-0490-438b-977b-b7d056973996", "2");
ma_hit_sub(min_dp, src, n_read, readLen, mini_overlap_length, cov);
// prt_dbg_rid_ovlp(src, -1, (char*)"c85c2e91-0490-438b-977b-b7d056973996", "3");
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "init_0", src, NULL, NULL, NULL);
detect_chimeric_reads(src, n_read, readLen, *cov, asm_opt.max_ov_diff_final*2.0, ul, UL_COV_THRES);
// prt_dbg_rid_ovlp(src, -1, (char*)"c85c2e91-0490-438b-977b-b7d056973996", "4");
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "init_1", src, NULL, *cov, NULL);
ma_hit_cut(src, n_read, readLen, mini_overlap_length, cov);
// prt_dbg_rid_ovlp(src, -1, (char*)"c85c2e91-0490-438b-977b-b7d056973996", "5");
ma_hit_flt(src, n_read, *cov, max_hang_length, mini_overlap_length);
// prt_dbg_rid_ovlp(src, -1, (char*)"c85c2e91-0490-438b-977b-b7d056973996", "6");
ma_hit_contained_advance(src, n_read, *cov, ruIndex, max_hang_length, mini_overlap_length);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "init_2", src, ruIndex, *cov, NULL);
// prt_dbg_rid_ovlp(src, -1, (char*)"7b70a587-f56c-48ac-adf8-b98a67365063_2", "7");
if(!ul) {
sg = ma_sg_gen(src, n_read, *cov, max_hang_length, mini_overlap_length);
if(asm_opt.prt_dbg_gfa) prt_dbg_gfa(sg, "raw", *cov, src, ruIndex, max_hang_length, mini_overlap_length);
if(asm_opt.prt_dbg_gfa) prt_dbg_gfa(sg, "raw", *cov, src, ruIndex, max_hang_length, mini_overlap_length);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "init_3", NULL, ruIndex, NULL, sg);
asg_arc_del_trans(sg, gap_fuzz);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "init_4", NULL, ruIndex, NULL, sg);
} else {
ug_opt_t uopt;
sg = ma_sg_gen_ul(src, n_read, *cov, ruIndex, max_hang_length, mini_overlap_length, UL_COV_THRES);
@@ -39871,6 +39950,17 @@ 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, uint8_t *cmk, int debug_g)
{
// const char *dbg_ids[] = {
// "6152e53b-a2ec-4f23-898f-462c9dee8e0f",
// "c5bbf86e-3146-4757-9882-55a5baebf770",
// };
// const uint64_t dbg_ids_n[] = {
// 3177240,
// 2198541,
// };
// hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "test_0", sources, NULL, NULL, NULL);
char *o_file = get_outfile_name(output_file_name);
ma_sub_t *coverage_cut = *coverage_cut_ptr;
asg_t *sg = *sg_ptr;
@@ -39894,8 +39984,10 @@ ma_sub_t **coverage_cut_ptr, uint8_t *cmk, int debug_g)
///normalize_ma_hit_t_single_side(sources, n_read);
normalize_ma_hit_t_single_side_advance(sources, n_read, asm_opt.is_ont, cmk);
// hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "test_1", sources, NULL, NULL, NULL);
// normalize_ma_hit_t_single_side_advance_mult(sources, n_read, asm_opt.thread_num);
normalize_ma_hit_t_single_side_advance(reverse_sources, n_read, 0, cmk);
// hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "test_2", sources, NULL, NULL, NULL);
// normalize_ma_hit_t_single_side_advance_mult(reverse_sources, n_read, asm_opt.thread_num);
if (ha_opt_triobin(&asm_opt))
@@ -39925,6 +40017,21 @@ ma_sub_t **coverage_cut_ptr, uint8_t *cmk, int debug_g)
// prt_dbg_rid_ovlp(sources, 27087, NULL, "0-b");
// prt_specific_overlap(sources, 22233, 22235, "0-c");
// prt_specific_overlap(sources, 22235, 22233, "0-c");
///@brief debug
if (asm_opt.flag & HA_F_VERBOSE_GFA) {
write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF, cmk);
debug_gfa:;
}
// hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ini_s", sources, NULL, NULL, NULL);
// hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ini_r", reverse_sources, NULL, NULL, NULL);
// exit(1);
sg = gen_init_sg(min_dp, n_read, mini_overlap_length, max_hang_length, gap_fuzz, sources, readLen, ruIndex,
&b_mask_t, &coverage_cut, asm_opt.ar?&UL_INF:NULL, te);
// if(asm_opt.ar) exit(1);
@@ -39955,24 +40062,29 @@ ma_sub_t **coverage_cut_ptr, uint8_t *cmk, int debug_g)
free(unlean_name);
}
// if (asm_opt.flag & HA_F_VERBOSE_GFA) {
// write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF);
// debug_gfa:;
// }
/**
///@brief debug
if (asm_opt.flag & HA_F_VERBOSE_GFA) {
write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF, cmk);
debug_gfa:;
}
**/
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_z", sources, ruIndex, coverage_cut, sg);
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, te);
ul_clean_gfa(&uopt, sg, sources, reverse_sources, ruIndex, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
0.6, asm_opt.max_short_tip, gap_fuzz, &b_mask_t, !!asm_opt.ar, ha_opt_triobin(&asm_opt), UL_COV_THRES, cmk, o_file);
/**
///@brief debug
if (asm_opt.flag & HA_F_VERBOSE_GFA) {
write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF);
write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF, cmk);
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, te);
// 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,
@@ -40229,7 +40341,7 @@ long long bubble_dist, int read_graph, int write)
min_thres = asm_opt.max_short_tip + 1;
if (asm_opt.flag & HA_F_VERBOSE_GFA)
{
if(load_debug_graph(/**NULL**/&sg, &sources, /**NULL**/&coverage_cut, output_file_name, &reverse_sources, &ruIndex, &UL_INF))
if(load_debug_graph(NULL/**&sg**/, &sources, NULL/**&coverage_cut**/, output_file_name, &reverse_sources, &ruIndex, &UL_INF, (asm_opt.is_ont)?(&cmk):(NULL)))
{
fprintf(stderr, "debug gfa has been loaded\n");
@@ -40249,20 +40361,39 @@ long long bubble_dist, int read_graph, int write)
}
///debug_info_of_specfic_read("m64011_190830_220126/31720629/ccs", sources, reverse_sources, -1, "beg");
// const char *dbg_ids[] = {
// "6152e53b-a2ec-4f23-898f-462c9dee8e0f",
// "c5bbf86e-3146-4757-9882-55a5baebf770",
// };
// const uint64_t dbg_ids_n[] = {
// 3177240,
// 2198541,
// };
if (!(asm_opt.flag & HA_F_BAN_ASSEMBLY))
{
// if(asm_opt.is_ont) handle_chemical_r(asm_opt.thread_num, R_INF.total_reads);
if(asm_opt.is_ont) {
uint64_t i = 0, j = 0;
// uint64_t i = 0, j = 0;
memset(ruIndex.index, -1, sizeof(uint32_t)*(ruIndex.len));
for (i = 0; i < n_read; i++) {
for (j = 0; j < sources[i].length; j++) sources[i].buffer[j].del = 0;
for (j = 0; j < reverse_sources[i].length; j++) reverse_sources[i].buffer[j].del = 0;
}
inital_ovlap(sources, reverse_sources, n_read, asm_opt.thread_num);
// hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ttst_0", sources, NULL, NULL, NULL);
// hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "trev_0", reverse_sources, NULL, NULL, NULL);
if(asm_opt.is_ont) cmk = gen_chemical_arc_rf(asm_opt.thread_num, R_INF.total_reads);
// hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ttst_1", sources, NULL, NULL, NULL);
if(asm_opt.del_hf) clean_arc_rf(asm_opt.thread_num, R_INF.total_reads);
// hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ttst_2", sources, NULL, NULL, NULL);
}
try_rescue_overlaps(sources, reverse_sources, n_read, 4, asm_opt.is_ont);
// hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ttst_3", sources, NULL, NULL, NULL);
clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length,
max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
+629 -62
View File
File diff suppressed because it is too large Load Diff
+176 -2
View File
@@ -411,6 +411,153 @@ void init_integer_ml_t(integer_ml_t *x, ul_resolve_t *u, uint64_t n_thread)
int32_t if_sup_chimeric(ma_hit_t_alloc* src, uint64_t rLen, asg64_v *b, int if_exact);
typedef struct {
const char *nn;
uint64_t nid;
} dbg_prt_t;
static void qry_nn_id(void *data, long i, int tid) // callback for kt_for()
{
dbg_prt_t *sl = (dbg_prt_t *)data; uint64_t ri = sl[i].nid, rl = 0; const char *nn = sl[i].nn;
rl = strlen(nn);
if((ri < R_INF.total_reads) && (Get_NAME_LENGTH((R_INF), ri) == rl) && (memcmp(nn, Get_NAME((R_INF), ri), Get_NAME_LENGTH((R_INF), ri)) == 0)) return;
for (ri = 0; ri < R_INF.total_reads; ri++) {
if((Get_NAME_LENGTH((R_INF), ri) == rl) && (memcmp(nn, Get_NAME((R_INF), ri), Get_NAME_LENGTH((R_INF), ri)) == 0)) {
sl[i].nid = ri;
return;
}
}
sl[i].nid = (uint64_t)-1;
}
void hc_dbg_prt_ma_hit_t(const char *ids0[], uint64_t n0, const uint64_t iids0[], uint64_t ni0, const char* cmd, ma_hit_t_alloc *paf, R_to_U* ruIndex, ma_sub_t* cov, asg_t* g)
{
if((!paf) && (!ruIndex) && (!cov) && (!g)) return;
uint64_t k, z, n, ni, rid, uri, n_thre = asm_opt.thread_num; dbg_prt_t *ri = NULL; const char **ids = NULL; const uint64_t *iids = NULL;
if (cmd == NULL) cmd = "";
fprintf(stderr, "\n[M::%s::%s]\t********************************\n", __func__, cmd);
const char *ids1[] = {
///"m84039_230928_213653_s3/181277190/ccs",///tip 0, contained read issue
///"7e61bf83-7f09-49df-88b9-2a864587b19f", ///tip1 -> 916fced8-524e-4a65-b04d-5b5ea6937e6e(qn::199123)
///phasing error -> ONT specific bias && HPC intoduce fake SNPs
///"43691917-c944-4102-8a01-d63a9c635dc3",///tip3 -> 29ee66bc-ef49-42ae-b7b4-15c37d46beac(qn::3373007) overcorrected in the first round
//when only 2 reads supported (the read itself has high sequencing error rate so hard to be aligned)
//we may also consider increase the error rate threshold
///"21b83ce9-9b94-445d-93b0-b9e83aaa4326",///tip3
"6f240062-2451-4212-a22d-abecc73ae0d9",///tip4->63a321ea-1103-4ca5-97e1-531e55f6bc84(qn::3646295)
///phasing error -> ONT specific bias && HPC intoduce fake SNPs && the DP was to nice for the following SNPs, which we need to fix as high-sequencing depth we have more such cases
// +[M::gen_rphase_dp0_single_path_hybrid_0_multi] rn::2 hf_only::0
// +[M::gen_rphase_dp0_single_path_hybrid_0_multi] site::175009 sc::2 n0::57 n1::3 rn::2 krn::2 b0l::16 b0h::41 b1l::0 b1h::3
// +[M::gen_rphase_dp0_single_path_hybrid_0_multi] site::82558 sc::-1 n0::80 n1::4 rn::2 krn::2 b0l::4 b0h::76 b1l::2 b1h::2
"a72d3f7f-d74f-48ab-a3f7-5ecc74ad08ff",///tip5, contained read issue
};
if(ids0) {
ids = ids0; n = n0;
} else {
ids = ids1; n = sizeof(ids1) / sizeof(ids1[0]);
}
const uint64_t iids1[] = {
///5038496,///tip 0, contained read issue
///3799659,///tip 1, contained read issue
///3313150,///tip 2, contained read issue
///4192664,///tip3, contained read issue
508213,///tip4
3203512,///tip5
};
if(iids0) {
iids = iids0; ni = ni0;
} else {
iids = iids1; ni = sizeof(iids1) / sizeof(iids1[0]);
}
CALLOC(ri, n);
for (k = uri = 0; k < n; k++) {
rid = ((k<ni)?(iids[k]):((uint64_t)-1));
if(rid < R_INF.total_reads) {
if ((Get_NAME_LENGTH((R_INF), rid) != strlen(ids[k])) || (memcmp(ids[k], Get_NAME((R_INF), rid), Get_NAME_LENGTH((R_INF), rid)))) {
rid = (uint64_t)-1;
}
}
if(rid >= R_INF.total_reads) uri++;
ri[k].nid = rid; ri[k].nn = ids[k];
}
if(uri) {
if(n_thre > n) n_thre = n;
kt_for(n_thre, qry_nn_id, ri, n);
}
ma_hit_t *h = NULL; uint64_t nv, v, w; asg_arc_t *av; uint32_t cid, is_u;
for (k = 0; k < n; k++) {
if(ri[k].nid != ((uint64_t)-1)) {
fprintf(stderr, "\n[M::%s::%s]\trid::%lu\t%.*s\trlen::%lu\n", __func__, cmd, ri[k].nid, (int)Get_NAME_LENGTH(R_INF, ri[k].nid), Get_NAME(R_INF, ri[k].nid), Get_READ_LENGTH(R_INF, ri[k].nid));
if(paf) {
fprintf(stderr, "[M::%s::%s]\tpaf::(abnm->%u\tfcor->%u)\n", __func__, cmd, paf[ri[k].nid].is_abnormal, paf[ri[k].nid].is_fully_corrected);
for (z = 0; z < paf[ri[k].nid].length; z++) {
h = &(paf[ri[k].nid].buffer[z]);
fprintf(stderr, "[%s::paf]\t%.*s(qn::%u)\t%u\t%u\t%u\t%c\t%.*s(tn::%u)\t%u\t%u\t%u\tml::%u\tbl::%u\tel::%u\tdel::%u\tnl_idl::%u", cmd, (int)Get_NAME_LENGTH(R_INF, Get_qn(*h)), Get_NAME((R_INF), Get_qn(*h)), Get_qn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_qn(*h)), Get_qs(*h), Get_qe(*h), "+-"[h->rev],
(int)Get_NAME_LENGTH(R_INF, Get_tn(*h)), Get_NAME((R_INF), Get_tn(*h)), Get_tn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_tn(*h)), Get_ts(*h), Get_te(*h), h->ml, h->bl,
h->el, h->del, h->no_l_indel);
if(ruIndex) {
get_R_to_U(ruIndex, Get_tn(*h), &cid, &is_u);
if((is_u != (uint32_t)-1) && (is_u == 0)) {
fprintf(stderr, "\tis_contain::1\tis_cr::1(crid->::%u::%.*s)", cid, (int)Get_NAME_LENGTH(R_INF, cid), Get_NAME(R_INF, cid));
}
}
fprintf(stderr, "\n");
}
}
if(g) {
fprintf(stderr, "[M::%s::%s]\tsg::(del->%u\tc->%u)\n", __func__, cmd, g->seq[ri[k].nid].del, g->seq[ri[k].nid].c);
v = ri[k].nid<<1;
nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
for (z = 0; z < nv; ++z) {
w = av[z].v;
fprintf(stderr, "[M::%s::sg]\t%.*s(%c)\t%.*s(%c)\tol::%u\tel::%u\tdel::%u\tstrong::%u\tnl_idl::%u\n", __func__,
(int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), "+-"[v&1],
(int)Get_NAME_LENGTH(R_INF, (w>>1)), Get_NAME(R_INF, (w>>1)), "+-"[w&1], av[z].ol, av[z].el, av[z].del, av[z].strong, av[z].no_l_indel);
}
v = (ri[k].nid<<1) + 1;
nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
for (z = 0; z < nv; ++z) {
w = av[z].v;
fprintf(stderr, "[M::%s::sg]\t%.*s(%c)\t%.*s(%c)\tol::%u\tel::%u\tdel::%u\tstrong::%u\tnl_idl::%u\n", __func__,
(int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), "+-"[v&1],
(int)Get_NAME_LENGTH(R_INF, (w>>1)), Get_NAME(R_INF, (w>>1)), "+-"[w&1], av[z].ol, av[z].el, av[z].del, av[z].strong, av[z].no_l_indel);
}
}
if(cov) {
fprintf(stderr, "[M::%s::%s]\tcov::(del->%u\tc->%u)\n", __func__, cmd, cov[ri[k].nid].del, cov[ri[k].nid].c);
}
if(ruIndex) {
fprintf(stderr, "[M::%s::%s]\t", __func__, cmd);
get_R_to_U(ruIndex, ri[k].nid, &cid, &is_u);
if((is_u != (uint32_t)-1) && (is_u == 0)) {
fprintf(stderr, "is_contain::1\tis_cr::1(rid->::%u::%.*s)\tis_cu::0\n", cid, (int)Get_NAME_LENGTH(R_INF, cid), Get_NAME(R_INF, cid));
} else if(is_u != (uint32_t)-1) {
fprintf(stderr, "is_contain::1\tis_cr::0\tis_cu::1(uid->::%u)\n", cid);
} else {
fprintf(stderr, "is_contain::0\tis_cr::0\tis_cu::0\n");
}
}
} else {
fprintf(stderr, "\n[M::%s::%s]\trid::-1\t%.*s\n", __func__, cmd, (int)strlen(ri[k].nn), ri[k].nn);
}
}
free(ri);
}
void print_edge(asg_arc_t *t, const char *cmd)
{
uint32_t v = t->ul>>32, w = t->v;
@@ -3057,41 +3204,53 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
if(asm_opt.is_ont) asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_a", NULL, rI, NULL, sg);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);///p_telo
// fprintf(stderr, "[M::%s] count_edges_v_w(sg, 49778, 49847)->%ld\n", __func__, count_edges_v_w(sg, 49778, 49847));
// if(is_ou) dedup_contain_g(uopt, sg);
// debug_info_of_specfic_node("c7ecbd6b-e09d-4042-93ac-2400839feaf6", sg, rI, "beg-1");
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_b", NULL, rI, NULL, sg);
for (i = 0; i < clean_round; i++, drop += step) {
if(drop > max_ovlp_drop_ratio) drop = max_ovlp_drop_ratio;
if(is_ou) {
if(drop <= 0.500001) min_diff = step_diff>>1;
else min_diff = step_diff;
}
// fprintf(stderr, "\n(0):i->%ld, drop->%f\n", i, drop);
if(asm_opt.is_ont) {
asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_0", NULL, rI, NULL, sg);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_1", NULL, rI, NULL, sg);
}
// fprintf(stderr, "(0):i->%ld, drop->%f\n", i, drop);
// prt_specfic_sge(sg, 22708, 22646, "--0--");
// print_vw_edge(sg, 34156, 34090, "0");
// stats_chimeric(sg, src, &bu);
if(!is_ou) asg_iterative_semi_circ(sg, src, &bu, max_tip, 1, uopt->te);///p_telo
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_3", NULL, rI, NULL, sg);
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1);
asg_arc_cut_chimeric(sg, src, &bu, is_ou?ou_thres:(uint32_t)-1, uopt->te);///p_telo
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_4", NULL, rI, NULL, sg);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_5", NULL, rI, NULL, sg);
// prt_specfic_sge(sg, 22708, 22646, "--1--");
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0);
asg_arc_cut_inexact(sg, src, &bu, max_tip, is_ou, is_trio, min_diff, ou_drop_rate/**, NULL**//**&dbg**/);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_6", NULL, rI, NULL, sg);
// debug_edges(&dbg, d, 2);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_7", NULL, rI, NULL, sg);
// prt_specfic_sge(sg, 22708, 22646, "--2--");
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1);
asg_arc_cut_length(sg, &bu, max_tip, drop, ou_drop_rate, is_ou, is_trio, 1, min_diff, 1, NULL, NULL, NULL);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_8", NULL, rI, NULL, sg);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_9", NULL, rI, NULL, sg);
// if(i == 0) {
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty3.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
@@ -3107,11 +3266,14 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1);
asg_arc_cut_bub_links(sg, &bu, HARD_OL_DROP, HARD_OL_SEC_DROP, HARD_OU_DROP, is_ou, asm_opt.large_pop_bubble_size, rev, rI, max_tip);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_10", NULL, rI, NULL, sg);
// prt_specfic_sge(sg, 22708, 22646, "--4--");
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1);
asg_arc_cut_complex_bub_links(sg, &bu, HARD_OL_DROP, HARD_OU_DROP, is_ou, b_mask_t);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_11", NULL, rI, NULL, sg);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_12", NULL, rI, NULL, sg);
// prt_specfic_sge(sg, 22708, 22646, "--5--");
/**
@@ -3133,7 +3295,9 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
if(asm_opt.is_ont) {
asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_0", NULL, rI, NULL, sg);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_1", NULL, rI, NULL, sg);
}
if(is_ou) min_diff = step_diff;
@@ -3147,6 +3311,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
if(clean_contain_g(uopt, sg, 1)) update_sg_uo(sg, src);
}
if(!is_ou) asg_iterative_semi_circ(sg, src, &bu, max_tip, 1, uopt->te);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_2", NULL, rI, NULL, sg);
// prt_specfic_sge(sg, 22708, 22646, "--sb-0---");
@@ -3155,23 +3320,29 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0);
asg_cut_large_indel(sg, &bu, max_tip, HARD_OU_DROP, is_ou, min_diff);///shoule we ignore ou here?
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_3", NULL, rI, NULL, sg);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_4", NULL, rI, NULL, sg);
// prt_specfic_sge(sg, 22708, 22646, "--sb-1---");
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty6.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "hybrid.dirty6.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
if(!is_ou) {
///asg_arc_del_triangular_directly might be unnecessary
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0);
asg_arc_cut_length(sg, &bu, max_tip, HARD_ORTHOLOGY_DROP/**min_ovlp_drop_ratio**/, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, 1, rev, rI, NULL);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_5", NULL, rI, NULL, sg);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_6", NULL, rI, NULL, sg);
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty7.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0);
asg_arc_cut_length(sg, &bu, max_tip, min_ovlp_drop_ratio, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, 1, rev, rI, &l_drop);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_7", NULL, rI, NULL, sg);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_8", NULL, rI, NULL, sg);
} else {
min_diff = step_diff; l_drop = 6000;
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0);
@@ -3184,6 +3355,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty8.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
if(!is_ou) asg_cut_semi_circ(sg, LIM_LEN, 1);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_9", NULL, rI, NULL, sg);
// prt_specfic_sge(sg, 22708, 22646, "--sb-3---");
/**
@@ -3198,10 +3370,12 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
rescue_bubble_by_chain(sg, uopt->coverage_cut, src, rev, (asm_opt.max_short_tip*2), 0.15, 3, rI, 0.05, 0.9, uopt->max_hang, uopt->min_ovlp, 10, uopt->gap_fuzz, b_mask_t);
**/
post_rescue(uopt, sg, src, rev, rI, b_mask_t, is_ou, cmk);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_10", NULL, rI, NULL, sg);
// prt_specfic_sge(sg, 22708, 22646, "--sb-4---");
ug_ext_gfa(uopt, sg, ug_ext_len);
// hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_11", NULL, rI, NULL, sg);
// if(is_ou) dedup_contain_g(uopt, sg);
+1
View File
@@ -41,5 +41,6 @@ void ug_ext_gfa(ug_opt_t *uopt, asg_t *sg, uint32_t max_len);
void update_sg_uo(asg_t *g, ma_hit_t_alloc *src);
uint32_t get_arcs(asg_t *g, uint32_t v, uint32_t* idx, uint32_t idx_n);
uint64_t ug_occ_w(uint64_t is, uint64_t ie, ma_utg_t *u);
void hc_dbg_prt_ma_hit_t(const char *ids0[], uint64_t n0, const uint64_t iids0[], uint64_t ni0, const char* cmd, ma_hit_t_alloc *paf, R_to_U* ruIndex, ma_sub_t* cov, asg_t* g);
#endif
+41 -7
View File
@@ -1420,7 +1420,7 @@ int load_ct_index(void **i_ct_idx, char* file_name)
return 1;
}
int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name, uint8_t force_rpaf_load)
{
char* gfa_name = (char*)malloc(strlen(file_name)+64);
if(r) sprintf(gfa_name, "%s.pt_flt", file_name);
@@ -1479,6 +1479,19 @@ int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t*
fwrite(&(r->paf[k].length), sizeof(r->paf[k].length), 1, fp);
fwrite(r->paf[k].buffer, sizeof((*(r->paf[k].buffer))), r->paf[k].length, fp);
}
int8_t rff = opt->post_syn;
if(force_rpaf_load) rff = 1;
fwrite(&rff, sizeof(rff), 1, fp);
if(rff) {
for (k = 0; k < r->total_reads; k++) {
fwrite(&(r->reverse_paf[k].is_fully_corrected), sizeof(r->reverse_paf[k].is_fully_corrected), 1, fp);
fwrite(&(r->reverse_paf[k].is_abnormal), sizeof(r->reverse_paf[k].is_abnormal), 1, fp);
fwrite(&(r->reverse_paf[k].length), sizeof(r->reverse_paf[k].length), 1, fp);
fwrite(r->reverse_paf[k].buffer, sizeof((*(r->reverse_paf[k].buffer))), r->reverse_paf[k].length, fp);
}
}
}
fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
@@ -1487,7 +1500,7 @@ int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t*
return 1;
}
int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_opt_t* opt, char* file_name)
int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_opt_t* opt, char* file_name, uint8_t force_rpaf_load)
{
char* gfa_name = (char*)malloc(strlen(file_name)+64);
if(r) sprintf(gfa_name, "%s.pt_flt", file_name);
@@ -1604,6 +1617,27 @@ int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_op
r->paf[k].buffer = (ma_hit_t*)malloc(sizeof(ma_hit_t)*r->paf[k].length);
fread(r->paf[k].buffer, sizeof((*(r->paf[k].buffer))), r->paf[k].length, fp);
}
int8_t rff = 0;
f_flag += fread(&rff, sizeof(rff), 1, fp);
if(!force_rpaf_load) {
opt->post_syn = rff;
fprintf(stderr, "[M::%s] overwritten post_syn to::%u\n", __func__, (uint32_t)opt->post_syn);
}
if(rff) {
for (k = 0; k < r->total_reads; k++) {
f_flag += fread(&(r->reverse_paf[k].is_fully_corrected), sizeof(r->reverse_paf[k].is_fully_corrected), 1, fp);
f_flag += fread(&(r->reverse_paf[k].is_abnormal), sizeof(r->reverse_paf[k].is_abnormal), 1, fp);
f_flag += fread(&(r->reverse_paf[k].length), sizeof(r->reverse_paf[k].length), 1, fp);
r->reverse_paf[k].size = r->reverse_paf[k].length;
r->reverse_paf[k].buffer = NULL;
if(r->reverse_paf[k].length == 0) continue;
r->reverse_paf[k].buffer = (ma_hit_t*)malloc(sizeof(ma_hit_t)*r->reverse_paf[k].length);
fread(r->reverse_paf[k].buffer, sizeof((*(r->reverse_paf[k].buffer))), r->reverse_paf[k].length, fp);
}
}
}
fclose(fp);
@@ -1619,16 +1653,16 @@ void refresh_pt_idx(void **flt_tab, ha_pt_t **ha_idx, All_reads *r, hifiasm_opt_
sprintf(gfa_name, "%s.ad", file_name);
if(is_w) {
write_pt_index(*flt_tab, *ha_idx, NULL, opt, gfa_name);
write_pt_index(*flt_tab, *ha_idx, NULL, opt, gfa_name, 0);
} else {
load_pt_index(flt_tab, ha_idx, NULL, opt, gfa_name);
load_pt_index(flt_tab, ha_idx, NULL, opt, gfa_name, 0);
}
free(gfa_name);
}
uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rcc, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load)
uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rcc, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load, uint8_t force_rpaf_load)
{
char* gfa_name = (char*)malloc(strlen(file_name)+64);
FILE *fp = NULL; int f_flag = 0; uint64_t rr0 = (uint64_t)-1, tot_rr0 = (uint64_t)-1;
@@ -1649,7 +1683,7 @@ uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rc
}
sprintf(gfa_name, "%s.r%lu", file_name, rr);
if(!load_pt_index(r_flt_tab, r_ha_idx, r, opt, gfa_name)) {
if(!load_pt_index(r_flt_tab, r_ha_idx, r, opt, gfa_name, force_rpaf_load)) {
destory_All_reads(r);
ha_pt_destroy(*r_ha_idx); (*r_ha_idx) = NULL;
ha_ft_destroy(*r_flt_tab); (*r_flt_tab) = NULL;
@@ -1673,7 +1707,7 @@ uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rc
write_cc_v(rcc, gfa_name);
sprintf(gfa_name, "%s.r%lu", file_name, rr);
write_pt_index(*r_flt_tab, *r_ha_idx, r, opt, gfa_name);
write_pt_index(*r_flt_tab, *r_ha_idx, r, opt, gfa_name, force_rpaf_load);
sprintf(gfa_name, "%s.r%lu.ht.bin", file_name, rr);
fp = fopen(gfa_name, "w");
+3 -3
View File
@@ -86,15 +86,15 @@ const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);
const ha_idxposl_t *ha_ptl_get(const ha_pt_t *h, uint64_t hash, int *n);
const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash);
int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name, uint8_t force_rpaf_load);
int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name, uint8_t force_rpaf_load);
void refresh_pt_idx(void **flt_tab, ha_pt_t **ha_idx, All_reads *r, hifiasm_opt_t *opt, char *file_name, uint8_t is_w);
int uidx_write(void *flt_tab, ha_pt_t *ha_idx, char* file_name, ma_ug_t *ug);
int uidx_load(void **r_flt_tab, ha_pt_t **r_ha_idx, char* file_name, ma_ug_t *ug);
int write_ct_index(void *ct_idx, char* file_name);
int load_ct_index(void **ct_idx, char* file_name);
int query_ct_index(void* ct_idx, uint64_t hash);
uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rcc, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load);
uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rcc, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load, uint8_t force_rpaf_load);
ha_abuf_t *ha_abuf_init_buf(void *km);
ha_abufl_t *ha_abufl_init_buf(void *km);
+8
View File
@@ -87,6 +87,14 @@ int main(int argc, char *argv[])
fprintf(stderr, "[M::%s::final] using %s mode\n", __func__, (asm_opt.simd_mm == 2) ? "AVX-512" : ((asm_opt.simd_mm == 1) ? "AVX2" : "non-SIMD"));
fprintf(stderr, "[M::%s::]\traw_aln::%d\tpost_syn::%d\n", __func__, asm_opt.realn_raw, asm_opt.post_syn);
if(asm_opt.recurrent_err_test) {
fprintf(stderr, "[M::%s::] Enable recurrent sequencing-error filtering\n", __func__);
} else {
fprintf(stderr, "[M::%s::] Disable recurrent sequencing-error filtering\n", __func__);
}
if(asm_opt.sec_in) ret = ha_assemble_pair();
else if(asm_opt.dbg_ovec_cal) ret = ha_ec_dbg();