mirror of
https://github.com/chhylp123/hifiasm.git
synced 2026-10-11 15:40:55 +08:00
Compare commits
2
Commits
f5078f7b23
...
bca2394b69
| Author | SHA1 | Date | |
|---|---|---|---|
 chhylp123
|
bca2394b69 | ||
|
chhylp123
|
7884b5ad88 |
+29
-14
@@ -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, &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,8 +2079,8 @@ 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;
|
||||
if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
|
||||
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());
|
||||
if (asm_opt.extract_list) {
|
||||
@@ -2093,16 +2096,26 @@ 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)) {
|
||||
for (r = asm_opt.number_of_round - 1; r >= 0; --r) {
|
||||
if(tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, 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, ((r==asm_opt.number_of_round)?(1):(0)))) {
|
||||
load_ct_index(&ha_ct_table, asm_opt.output_file_name); r0 = r;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if(r < 0) r = 0;
|
||||
if((asm_opt.dbg_ec_rr >= 0) && (asm_opt.dbg_ec_rr != r)) {
|
||||
fprintf(stderr, "[E::%s] no matching debug error-correction bins found\n", __func__);
|
||||
exit(1);
|
||||
}
|
||||
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
|
||||
@@ -2123,6 +2136,7 @@ int ha_assemble(void)
|
||||
// fprintf(stderr, "[M::%s] # bases: %lld; # corrected bases: %lld; # recorrected bases: %lld\n", __func__,
|
||||
// asm_opt.num_bases, asm_opt.num_corrected_bases, asm_opt.num_recorrected_bases);
|
||||
// fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0);
|
||||
if(asm_opt.dbg_ec_rr >= 0) exit(1);
|
||||
}
|
||||
if (asm_opt.flag & HA_F_WRITE_EC) {
|
||||
if(asm_opt.is_sc) Output_corrected_fastq();
|
||||
@@ -2131,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)) {
|
||||
@@ -2149,7 +2164,7 @@ int ha_assemble(void)
|
||||
build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf,
|
||||
R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round,
|
||||
asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded);
|
||||
destory_All_reads(&R_INF); if(asm_opt.dbg_bam) destroy_cc_v(&scb);
|
||||
destory_All_reads(&R_INF); /**if(asm_opt.dbg_bam)**/ destroy_cc_v(&scb);
|
||||
return 0;
|
||||
}
|
||||
|
||||
@@ -2186,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) {
|
||||
|
||||
+131
-2
@@ -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,6 +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},///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 }
|
||||
};
|
||||
@@ -105,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");
|
||||
@@ -139,6 +153,33 @@ 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);
|
||||
@@ -364,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;
|
||||
@@ -444,6 +485,24 @@ void init_opt(hifiasm_opt_t* asm_opt)
|
||||
asm_opt->simd_mm = -1;
|
||||
|
||||
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)
|
||||
@@ -834,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;
|
||||
}
|
||||
|
||||
@@ -1121,6 +1200,26 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
|
||||
asm_opt->simd_mm = atoi(opt.arg);
|
||||
} else if (c == 378) {
|
||||
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);
|
||||
}
|
||||
@@ -1162,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);
|
||||
}
|
||||
|
||||
+21
-1
@@ -5,7 +5,7 @@
|
||||
#include <pthread.h>
|
||||
#include <stdint.h>
|
||||
|
||||
#define HA_VERSION "0.25.1-r910"
|
||||
#define HA_VERSION "0.25.1-r933"
|
||||
|
||||
#define VERBOSE 0
|
||||
|
||||
@@ -204,6 +204,26 @@ typedef struct {
|
||||
|
||||
int8_t simd_mm;
|
||||
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;
|
||||
|
||||
+2465
-117
File diff suppressed because it is too large
Load Diff
@@ -1456,8 +1456,9 @@ void gen_hc_r_alin_adv(gen_hc_aln_t *ez);
|
||||
uint64_t gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, 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, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *kp, uint8_t *hpf);
|
||||
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);
|
||||
void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, 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);
|
||||
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 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);
|
||||
|
||||
+481
@@ -1540,6 +1540,63 @@ inline int32_t comput_sc_ch_ec(const k_mer_hit *ai, const k_mer_hit *aj, double
|
||||
return sc;
|
||||
}
|
||||
|
||||
inline int32_t comput_sc_ch_ec_global(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol, uint8_t no_adj)
|
||||
{
|
||||
///ai is the suffix of aj
|
||||
int32_t dq, dr, dd, dg, q_span, sc; double dg_of;
|
||||
dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
|
||||
dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
|
||||
if((dq < 0) || (dr < 0)) return INT32_MIN;
|
||||
if((no_adj) && ((dq == 0) || (dr == 0))) return INT32_MIN;
|
||||
|
||||
dd = dr > dq? dr - dq : dq - dr;//gap
|
||||
if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, sl, ol))) return INT32_MIN;
|
||||
|
||||
dg = dr < dq? dr : dq;//len
|
||||
q_span = ai->cnt&(0xffu);
|
||||
sc = q_span < dg? q_span : dg;
|
||||
sc = normal_w(sc, ((int32_t)(ai->cnt>>8)));
|
||||
if (dd || (dg > q_span && dg > 0)) {
|
||||
double lin_pen, a_pen;
|
||||
lin_pen = (chn_pen_gap*(double)dd);
|
||||
dg_of = (dg>0)?((double)dg):(0.333333);
|
||||
a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate);
|
||||
///for long gap
|
||||
// if(lin_pen > a_pen) lin_pen = a_pen;
|
||||
if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen));
|
||||
else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen));
|
||||
lin_pen += (chn_pen_skip*dg_of);
|
||||
sc -= (int32_t)lin_pen;
|
||||
}
|
||||
return sc;
|
||||
}
|
||||
|
||||
inline int32_t comput_sc_ff_adv(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol)
|
||||
{
|
||||
///ai is the suffix of aj
|
||||
int32_t dq, dr, dd, dg, q_span, sc;
|
||||
dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
|
||||
if(dq < 0) return INT32_MIN;
|
||||
dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
|
||||
if(dr < 0) return INT32_MIN;
|
||||
dd = dr > dq? dr - dq : dq - dr;//gap
|
||||
// if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, sl, ol))) return INT32_MIN;
|
||||
dg = dr < dq? dr : dq;//len
|
||||
if(dg <= 0) return INT32_MIN;
|
||||
q_span = ai->cnt&(0xffu);
|
||||
sc = q_span < dg? q_span : dg;
|
||||
sc = normal_w(sc, ((int32_t)(ai->cnt>>8)));
|
||||
if (dd || (dg > q_span && dg > 0)) {
|
||||
double lin_pen, a_pen;
|
||||
lin_pen = (chn_pen_gap*(double)dd);
|
||||
a_pen = ((double)(sc))*((((double)dd)/((double)dg))/bw_rate);
|
||||
if(lin_pen > a_pen) lin_pen = a_pen;
|
||||
lin_pen += (chn_pen_skip*(double)dg);
|
||||
sc -= (int32_t)lin_pen;
|
||||
}
|
||||
return sc;
|
||||
}
|
||||
|
||||
inline int32_t comput_sc_ff(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol)
|
||||
{
|
||||
///ai is the suffix of aj
|
||||
@@ -2283,6 +2340,315 @@ uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n,
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
void quick_ck_lchain_global(k_mer_hit* a, int64_t a_n, int64_t xl, int64_t yl, double chn_pen_gap, double chn_pen_skip, double bw_rate,
|
||||
int64_t *p, int64_t *t, int32_t *f, int32_t *ii, int32_t *gf, int64_t *gp, int64_t *si, int64_t *ei)
|
||||
{
|
||||
*si = 0; *ei = a_n; gf[0] = gf[1] = gp[0] = gp[1] = INT32_MIN;
|
||||
if((a_n <= 0) || (xl <= 0) || (yl <= 0)) return;
|
||||
int64_t l, k, is_srt = 1, z; k_mer_hit *ai, *aj, ft, fz;
|
||||
int64_t dq, dr, dd, dg, q_span, sc, csc, ddt; uint8_t ff; double lin_pen, a_pen, dg_of;
|
||||
ft.cnt = fz.cnt = 0xFFFFFF00u; ft.readID = fz.readID = a[0].readID;
|
||||
ft.offset = ft.self_offset = 0; fz.offset = yl-1; fz.self_offset = xl-1;
|
||||
|
||||
for (k = 1, l = 0; k <= a_n; k++) {
|
||||
if(k == a_n || a[k].strand != a[l].strand) {
|
||||
t[k-1] = 0; ii[k-1] = 0;
|
||||
|
||||
if(is_srt) {
|
||||
ddt = 0; ff = 0; p[l] = f[l] = INT32_MIN;
|
||||
|
||||
ft.strand = a[l].strand;
|
||||
aj = &ft; z = l; ai = &a[z];
|
||||
dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
|
||||
dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
|
||||
dd = dr > dq? dr - dq : dq - dr;//gap
|
||||
if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, xl, yl))) ff = 1;
|
||||
if(!ff) {
|
||||
dg = dr < dq? dr : dq;//len
|
||||
q_span = ai->cnt&(0xffu);
|
||||
sc = q_span < dg? q_span : dg;
|
||||
sc = normal_w(sc, ((int32_t)(ai->cnt>>8)));
|
||||
if (dd || (dg > q_span && dg > 0)) {
|
||||
lin_pen = (chn_pen_gap*(double)dd);
|
||||
dg_of = ((dg>0)?(dg):(0.333333));
|
||||
a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate);
|
||||
///for long gap
|
||||
// if(lin_pen > a_pen) lin_pen = a_pen;
|
||||
if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen));
|
||||
else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen));
|
||||
lin_pen += (chn_pen_skip*dg_of);
|
||||
sc -= (int32_t)lin_pen;
|
||||
}
|
||||
|
||||
csc = a[z].cnt&(0xffu); if(sc < csc) ff = 1;
|
||||
if(!ff) {
|
||||
p[z] = -1; f[z] = sc; ddt += dd;
|
||||
}
|
||||
}
|
||||
|
||||
if(!ff) {
|
||||
for (z = l + 1; z < k; z++) {
|
||||
///roughly same to comput_sc_ch(&a[z], &a[z-1])
|
||||
ai = &a[z]; aj = &a[z-1];
|
||||
dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
|
||||
if(dq <= 0) break;
|
||||
dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
|
||||
if(dr <= 0) break;
|
||||
dd = dr > dq? dr - dq : dq - dr;//gap
|
||||
|
||||
if((dd > 16) && (dd > cal_bw(&(a[z]), &(a[z-1]), bw_rate, xl, yl))) break;
|
||||
dg = dr < dq? dr : dq;//len
|
||||
q_span = ai->cnt&(0xffu);
|
||||
sc = q_span < dg? q_span : dg;
|
||||
sc = normal_w(sc, ((int32_t)(ai->cnt>>8)));
|
||||
if (dd || (dg > q_span && dg > 0)) {
|
||||
lin_pen = (chn_pen_gap*(double)dd);
|
||||
dg_of = ((dg>0)?(dg):(0.333333));
|
||||
a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate);
|
||||
///for long gap
|
||||
// if(lin_pen > a_pen) lin_pen = a_pen;
|
||||
if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen));
|
||||
else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen));
|
||||
lin_pen += (chn_pen_skip*dg_of);
|
||||
sc -= (int32_t)lin_pen;
|
||||
}
|
||||
|
||||
sc += f[z-1]; csc = a[z].cnt&(0xffu); if(sc < csc) break;
|
||||
p[z] = z - 1; f[z] = sc; ddt += dd;
|
||||
}
|
||||
if(z < k) ff = 1;
|
||||
}
|
||||
|
||||
if(!ff) {
|
||||
fz.strand = a[l].strand;
|
||||
ai = &fz; aj = &(a[k-1]);
|
||||
dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset);
|
||||
dr = (int64_t)(ai->offset) - (int64_t)(aj->offset);
|
||||
dd = dr > dq? dr - dq : dq - dr;//gap
|
||||
if(((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, xl, yl)))||
|
||||
(((ddt + dd) > 16) && ((ddt + dd) > cal_bw(&fz, &ft, bw_rate, xl, yl)))) {
|
||||
ff = 1;
|
||||
} else {
|
||||
dg = dr < dq? dr : dq;//len
|
||||
sc = q_span = 0;
|
||||
if (dd || (dg > q_span && dg > 0)) {
|
||||
lin_pen = (chn_pen_gap*(double)dd);
|
||||
dg_of = ((dg>0)?(dg):(0.333333));
|
||||
a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate);
|
||||
///for long gap
|
||||
// if(lin_pen > a_pen) lin_pen = a_pen;
|
||||
if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen));
|
||||
else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen));
|
||||
lin_pen += (chn_pen_skip*dg_of);
|
||||
sc -= (int32_t)lin_pen;
|
||||
}
|
||||
sc += f[k-1]; ///csc = a[z].cnt&(0xffu); if(sc < csc) break;
|
||||
gp[fz.strand] = k-1; gf[fz.strand] = sc; ///ddt += dd;
|
||||
if((*ei) > k) {
|
||||
(*si) = k;
|
||||
} else {
|
||||
(*ei) = l;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
l = k; is_srt = 1;
|
||||
} else {
|
||||
if((a[k].self_offset <= a[k-1].self_offset) || (a[k].offset <= a[k-1].offset)) is_srt = 0;
|
||||
t[k-1] = 0; ii[k-1] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
uint64_t lchain_qdp_global_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t des_idx,
|
||||
Chain_Data* dp, overlap_region_alloc* res, int64_t max_skip, int64_t max_iter,
|
||||
int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate,
|
||||
uint32_t xid, int64_t xl, int64_t yl, int64_t quick_check, uint32_t apend_be,
|
||||
int64_t gen_cigar, int64_t khit_n)
|
||||
{
|
||||
if(a_n <= 0) return 0;
|
||||
int64_t *p, *t, *gp, max_f, n_skip, st, max_j, end_j, sc, max_ii, ovl, min_sc, ch_n, si, ei;
|
||||
int32_t *f, max, tmp, *ii, *gf; int64_t i, k, j, cL = 0; k_mer_hit* a; k_mer_hit* des; k_mer_hit *swap, ft; overlap_region *z;
|
||||
resize_Chain_Data(dp, a_n + 2, NULL); ch_n = 1;
|
||||
t = dp->tmp; f = dp->score; p = dp->pre; ii = dp->occ; gp = p + a_n; gf = f + a_n; gf[0] = gf[1] = gp[0] = gp[1] = INT32_MIN;
|
||||
|
||||
a = cl->list + a_idx; des = cl->list + des_idx;
|
||||
if(quick_check) {
|
||||
quick_ck_lchain_global(a, a_n, xl, yl, chn_pen_gap, chn_pen_skip, bw_rate, p, t, f, ii, gf, gp, &si, &ei);
|
||||
} else {
|
||||
si = 0; ei = a_n; memset(t, 0, (a_n*sizeof((*t))));
|
||||
}
|
||||
|
||||
ft.cnt = 0xFFFFFF00u; ft.readID = a[0].readID; ft.strand = 0; ft.offset = ft.self_offset = 0;
|
||||
for (i = st = si, max_ii = -1; i < ei; ++i) {
|
||||
///max_f = a[i].cnt&(0xffu);
|
||||
ft.strand = a[i].strand;
|
||||
|
||||
|
||||
n_skip = 0; max_j = end_j = -1;
|
||||
if ((i-st) > max_iter) st = i-max_iter;
|
||||
while (a[i].strand != a[st].strand) ++st;
|
||||
|
||||
for (j = i - 1; j >= st; --j) {
|
||||
sc = comput_sc_ch_ec_global(&a[i], &a[j], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl);
|
||||
if (sc == INT32_MIN) continue;
|
||||
sc += f[j];
|
||||
if (sc > max_f) {
|
||||
max_f = sc, max_j = j;
|
||||
if (n_skip > 0) --n_skip;
|
||||
} else if (t[j] == (int32_t)i) {
|
||||
if (++n_skip > max_skip)
|
||||
break;
|
||||
}
|
||||
if (p[j] >= 0) t[p[j]] = i;
|
||||
}
|
||||
end_j = j;
|
||||
|
||||
if ((max_ii<0) || (a[i].self_offset>a[max_ii].self_offset+max_dis) || (a[i].strand!=a[max_ii].strand)) {
|
||||
max = INT32_MIN; max_ii = -1;
|
||||
for (j=i-1; (j>=st) && (a[i].self_offset<=max_dis+a[j].self_offset)&&(a[i].strand==a[j].strand); --j) {
|
||||
if (max < f[j]) {
|
||||
max = f[j], max_ii = j;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if ((max_ii >= 0) && (max_ii < end_j) && (a[i].strand == a[max_ii].strand)) {///just have a try with a[i]<->a[max_ii]
|
||||
tmp = comput_sc_ch_ec(&a[i], &a[max_ii], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl);
|
||||
if (tmp != INT32_MIN && max_f < tmp + f[max_ii])
|
||||
max_f = tmp + f[max_ii], max_j = max_ii;
|
||||
}
|
||||
f[i] = max_f; p[i] = max_j;
|
||||
if ((max_ii < 0) || ((a[i].self_offset<=max_dis+a[max_ii].self_offset)&&(a[i].strand==a[max_ii].strand)&&(f[max_ii]<f[i]))) {
|
||||
max_ii = i;
|
||||
}
|
||||
if(f[i] >= msc) {
|
||||
ovl = get_chainLen(a[i].self_offset, a[i].self_offset, xl, a[i].offset, a[i].offset, yl);
|
||||
if(f[i] > msc || ovl < movl) {
|
||||
msc = f[i]; msc_i = i; movl = ovl;
|
||||
}
|
||||
}
|
||||
if(f[i] < plus) plus = f[i];
|
||||
ii[i] = 0;///for mcopy, not here
|
||||
// if(a_n && (a[0].readID == 27105 || a[0].readID == 7603)) {///r833
|
||||
// fprintf(stderr, "i::%ld[M::%s::rid->%u::%c] q::%u, t::%u, st::%ld, max_ii::%ld, f[i]::%d, p[i]::%ld, msc_i::%ld, msc::%ld, movl::%ld\n",
|
||||
// i, __func__, a[i].readID, "+-"[a[i].strand],
|
||||
// a[i].self_offset, a[i].offset, st, max_ii, f[i], p[i], msc_i, msc, movl);
|
||||
// }
|
||||
}
|
||||
|
||||
for (i = msc_i, cL = 0; i >= 0; i = p[i]) { ii[i] = 1; t[cL++] = i;}///label the best chain
|
||||
|
||||
if(mcopy_num > 1) {
|
||||
// if(a[0].readID == 4412344) {
|
||||
// fprintf(stderr, "[M::%s::] msc::%ld, cL::%ld\n", __func__, msc, cL);
|
||||
// }
|
||||
if(cL >= mcopy_khit_cutoff) {///if there are too few k-mers, disable mcopy
|
||||
msc -= plus; min_sc = msc*mcopy_rate; ii[msc_i] = 0;
|
||||
for (i = ch_n = 0; i < a_n; ++i) {///make all f[] positive
|
||||
f[i] -= plus; if(i >= ch_n) t[i] = 0;
|
||||
if((!(ii[i])) && (f[i] >= min_sc)) {///!(ii[i]): skip the best chain
|
||||
t[ch_n] = ((uint64_t)f[i])<<32; t[ch_n] += (i<<1); ch_n++;
|
||||
}
|
||||
}
|
||||
// if(a[0].readID == 4412344) {
|
||||
// fprintf(stderr, "[M::%s::] msc::%ld, min_sc::%ld, cL::%ld, ch_n::%ld, mcopy_num::%ld\n", __func__, msc, min_sc, cL, ch_n, mcopy_num);
|
||||
// }
|
||||
if(ch_n > 1) {
|
||||
int64_t n_v, n_v0, ni, n_u, n_u0 = res->length;
|
||||
radix_sort_hc64i(t, t + ch_n);
|
||||
for (k = ch_n-1, n_v = n_u = 0; k >= 0 && n_u < mcopy_num; --k) {
|
||||
n_v0 = n_v;
|
||||
for (i = ((uint32_t)t[k])>>1; i >= 0 && (t[i]&1) == 0; ) {
|
||||
ii[n_v++] = i; t[i] |= 1; i = p[i];
|
||||
}
|
||||
if(n_v0 == n_v) continue;
|
||||
sc = (i<0?(t[k]>>32):((t[k]>>32)-f[i]));
|
||||
// if(a[0].readID == 4412344) {
|
||||
// fprintf(stderr, "+[M::%s::] sc::%ld, n_a::%ld\n", __func__, sc, n_v-n_v0);
|
||||
// }
|
||||
if(sc >= min_sc) {
|
||||
kv_pushp_ol(overlap_region, (*res), &z);
|
||||
push_ovlp_chain_qgen(z, xid, xl, yl, sc+plus, &(a[ii[n_v-1]]), &(a[ii[n_v0]]));
|
||||
// if(a[0].readID == 4412344) {
|
||||
// fprintf(stderr, "-[M::%s::] sc::%ld, n_a::%ld, q::[%u,%u), t::[%u,%u), %c\n", __func__, sc, n_v-n_v0, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, "+-"[z->y_pos_strand]);
|
||||
// }
|
||||
///mcopy_khit_cutoff <= 1: disable the mcopy_khit_cutoff filtering, for the realignment
|
||||
// if((mcopy_khit_cutoff <= 1) || ((z->x_pos_e+1-z->x_pos_s) <= (movl<<2))) {
|
||||
if((!n_u) || (n_v - n_v0 > 1)) {
|
||||
z->align_length = n_v-n_v0; z->x_id = n_v0;
|
||||
n_u++;
|
||||
} else {///non-best is tiny
|
||||
res->length--; n_v = n_v0;
|
||||
}
|
||||
} else {
|
||||
n_v = n_v0;
|
||||
}
|
||||
}
|
||||
|
||||
// if(n_u > 1) ks_introsort_or_sss(n_u, res->list + n_u0);
|
||||
// res->length = n_u0 + filter_non_ovlp_xchains(res->list + n_u0, n_u, &n_v);
|
||||
n_u = res->length;
|
||||
if(n_u > n_u0 + 1) {
|
||||
kv_resize_cl(k_mer_hit, (*cl), (n_v+cl->length));
|
||||
a = cl->list + a_idx; des = cl->list + des_idx; swap = cl->list + cl->length;
|
||||
for (k = n_u0, i = n_v0 = n_v = 0; k < n_u; k++) {
|
||||
z = &(res->list[k]);
|
||||
z->non_homopolymer_errors = des_idx + i;
|
||||
n_v0 = z->x_id; ni = z->align_length;
|
||||
for (j = 0; j < ni; j++, i++) {
|
||||
///k0 + (ni - j - 1)
|
||||
swap[i] = a[ii[n_v0 + (ni- j - 1)]];
|
||||
swap[i].readID = k;
|
||||
}
|
||||
z->x_id = xid;
|
||||
if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, swap+i-ni, ni);
|
||||
if(!khit_n) z->align_length = 0;
|
||||
}
|
||||
memcpy(des, swap, i*sizeof((*swap))); //assert(i == ch_n);
|
||||
|
||||
// fprintf(stderr, "[M::%s::msc->%ld] msc_k_hits::%u, cL::%ld, min_sc::%ld, best_sc::%ld, n_u0_sc::%d, mcopy_rate::%f, # chains::%ld\n",
|
||||
// __func__, msc, res->list[n_u0].align_length, cL, min_sc, msc+plus, res->list[n_u0].shared_seed,
|
||||
// mcopy_rate, n_u-n_u0);
|
||||
} else if(n_u == n_u0 + 1) {
|
||||
z = &(res->list[n_u0]); k = n_u0; i = 0;
|
||||
z->non_homopolymer_errors = des_idx + i;
|
||||
n_v0 = z->x_id; ni = z->align_length;
|
||||
for (j = 0; j < ni; j++, i++) {
|
||||
///k0 + (ni - j - 1)
|
||||
des[i] = a[ii[n_v0 + (ni- j - 1)]];
|
||||
des[i].readID = k;
|
||||
}
|
||||
z->x_id = xid;
|
||||
if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, des+i-ni, ni);
|
||||
if(!khit_n) z->align_length = 0;
|
||||
}
|
||||
return i;
|
||||
} else {
|
||||
msc += plus; i = msc_i; cL = 0;
|
||||
while (i >= 0) {t[cL++] = i; i = p[i];}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
///a[] has been sorted by self_offset
|
||||
// i = msc_i; cL = 0;
|
||||
// while (i >= 0) {t[cL++] = i; i = p[i];}
|
||||
kv_pushp_ol(overlap_region, (*res), &z);
|
||||
push_ovlp_chain_qgen(z, xid, xl, yl, msc, &(a[t[cL-1]]), &(a[t[0]]));
|
||||
for (i = 0; i < cL; i++) {des[i] = a[t[cL-i-1]]; des[i].readID = res->length-1;}
|
||||
z->non_homopolymer_errors = des_idx;
|
||||
if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, des, cL);
|
||||
if(khit_n) z->align_length = cL;
|
||||
return cL;
|
||||
}
|
||||
**/
|
||||
|
||||
|
||||
#define rev_khit(an, xl, yl) do { \
|
||||
(an).self_offset = (xl)-1-((an).self_offset+1-((an).cnt&((uint32_t)(0xffu)))); \
|
||||
(an).offset = (yl)-1-((an).offset+1-((an).cnt&((uint32_t)(0xffu))));\
|
||||
@@ -2453,6 +2819,121 @@ uint64_t lchain_qdp_fix(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_s
|
||||
}
|
||||
|
||||
|
||||
uint64_t lchain_qdp_fix_adv(k_mer_hit *a, int64_t a_n, Chain_Data* dp, int64_t max_skip,
|
||||
int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip,
|
||||
double bw_rate, int64_t xl, int64_t yl, int64_t quick_check,
|
||||
int64_t left_fix, int64_t right_fix, k_mer_hit *res)
|
||||
{
|
||||
if(a_n <= 0) return 0;
|
||||
|
||||
int64_t *p, *t, max_f, n_skip, st, max_j, end_j, sc, msc, msc_i, bw, max_ii, ovl, movl;
|
||||
int32_t *f, max, tmp; int64_t i, j, ret, cL = 0;
|
||||
resize_Chain_Data(dp, a_n, NULL);
|
||||
t = dp->tmp; f = dp->score; p = dp->pre;
|
||||
bw = ((xl < yl)?xl:yl); bw *= bw_rate;
|
||||
msc = msc_i = -1; movl = INT32_MAX;
|
||||
|
||||
if(quick_check) {
|
||||
ret = lchain_qcheck(a, a_n, dp, bw_rate);
|
||||
if (ret > 0) {
|
||||
a_n = ret; msc_i = a_n-1; msc = f[msc_i];
|
||||
goto skip_ldp;
|
||||
}
|
||||
}
|
||||
|
||||
memset(t, 0, (a_n*sizeof((*t))));
|
||||
for (i = st = 0, max_ii = -1; i < a_n; ++i) {
|
||||
max_f = a[i].cnt&(0xffu); if(left_fix && i > 0) max_f = INT32_MIN;
|
||||
n_skip = 0; max_j = end_j = -1;
|
||||
if ((i-st) > max_iter) st = i-max_iter;
|
||||
|
||||
for (j = i - 1; j >= 0; --j) {
|
||||
if(left_fix && f[j] == INT32_MIN)continue;
|
||||
sc = comput_sc_ff_adv(&a[i], &a[j], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl);
|
||||
if (sc == INT32_MIN) continue;
|
||||
sc += f[j];
|
||||
if (sc > max_f) {
|
||||
max_f = sc, max_j = j;
|
||||
if (n_skip > 0) --n_skip;
|
||||
} else if (t[j] == (int32_t)i) {
|
||||
if ((++n_skip) > max_skip) {
|
||||
if((max_j != -1) || (left_fix == 0)) break;
|
||||
}
|
||||
}
|
||||
if (p[j] >= 0) t[p[j]] = i;
|
||||
///put it here will allow at least one prefix no matter max_dis
|
||||
///this is special for gap filling, not for chaining
|
||||
if (a[i].self_offset > (max_dis + a[j].self_offset)) {
|
||||
if((max_j != -1)) break;
|
||||
}
|
||||
if (j < st) {
|
||||
if((max_j != -1) || (left_fix == 0)) break;
|
||||
}
|
||||
}
|
||||
end_j = j;
|
||||
|
||||
if (max_ii < 0 || ((int64_t)a[i].self_offset) - ((int64_t)a[max_ii].self_offset) > max_dis) {
|
||||
max = INT32_MIN; max_ii = -1;
|
||||
for (j = i - 1; (j >= st) && ((((int64_t)a[i].self_offset)-((int64_t)a[j].self_offset))<=max_dis); --j) {
|
||||
if ((f[j] != INT32_MIN) && (max < f[j])) {
|
||||
max = f[j], max_ii = j;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if ((max_ii >= 0) && (max_ii < end_j) && (f[max_ii] != INT32_MIN)) {///just have a try with a[i]<->a[max_ii]
|
||||
tmp = comput_sc_ff_adv(&a[i], &a[max_ii], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl);
|
||||
if (tmp != INT32_MIN && max_f < tmp + f[max_ii])
|
||||
max_f = tmp + f[max_ii], max_j = max_ii;
|
||||
}
|
||||
|
||||
p[i] = max_j; f[i] = max_f;
|
||||
|
||||
if ((max_ii < 0) || (((((int64_t)a[i].self_offset)-((int64_t)a[max_ii].self_offset))<=max_dis) && (f[max_ii]<f[i]))) {
|
||||
max_ii = i;
|
||||
}
|
||||
|
||||
if(f[i] >= msc) {
|
||||
ovl = get_chainLen(a[i].self_offset, a[i].self_offset, xl, a[i].offset, a[i].offset, yl);
|
||||
if(f[i] > msc || ovl < movl) {
|
||||
msc = f[i]; msc_i = i; movl = ovl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
skip_ldp:
|
||||
if(right_fix && f[a_n-1] == INT32_MIN) return 0;
|
||||
if(right_fix) msc_i = a_n-1;
|
||||
///a[] has been sorted by self_offset
|
||||
i = msc_i; cL = 0;
|
||||
while (i >= 0) {
|
||||
t[cL++] = i; msc_i = i; i = p[i];
|
||||
}
|
||||
|
||||
n_skip = cL>>1;
|
||||
for (i = 0; i < n_skip; i++) {
|
||||
msc_i = t[i]; t[i] = t[cL-i-1]; t[cL-i-1] = msc_i;
|
||||
}
|
||||
|
||||
if((cL > 0) && (right_fix) && (t[cL-1] != (a_n-1))) {
|
||||
cL = 0;
|
||||
}
|
||||
|
||||
if((cL > 0) && (left_fix) && (t[0] != 0)) {
|
||||
cL = 0;
|
||||
}
|
||||
|
||||
if(cL > 0 && res) {
|
||||
for (i = 0; i < cL; i++) {
|
||||
res[i] = a[t[i]];
|
||||
}
|
||||
}
|
||||
|
||||
return cL;
|
||||
}
|
||||
|
||||
|
||||
|
||||
uint64_t lchain_refine(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp,
|
||||
int64_t max_skip, int64_t max_iter, int64_t max_dis, int64_t long_gap)
|
||||
{
|
||||
|
||||
@@ -234,6 +234,10 @@ uint64_t lchain_qdp_fix(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_s
|
||||
int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip,
|
||||
double bw_rate, int64_t xl, int64_t yl, int64_t quick_check,
|
||||
int64_t left_fix, int64_t right_fix);
|
||||
uint64_t lchain_qdp_fix_adv(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_skip,
|
||||
int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip,
|
||||
double bw_rate, int64_t xl, int64_t yl, int64_t quick_check,
|
||||
int64_t left_fix, int64_t right_fix, k_mer_hit *res);
|
||||
uint64_t lchain_simple(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp,
|
||||
int64_t max_skip, int64_t max_iter);
|
||||
uint64_t lchain_simple0(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, int64_t max_skip, int64_t max_iter);
|
||||
@@ -253,6 +257,12 @@ uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n,
|
||||
int64_t gen_cigar, int64_t enable_mcopy, double mcopy_rate, int64_t mcopy_khit_cutoff,
|
||||
int64_t khit_n);
|
||||
|
||||
uint64_t lchain_qdp_global_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t des_idx,
|
||||
Chain_Data* dp, overlap_region_alloc* res, int64_t max_skip, int64_t max_iter,
|
||||
int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate,
|
||||
uint32_t xid, int64_t xl, int64_t yl, int64_t quick_check, uint32_t apend_be,
|
||||
int64_t gen_cigar, int64_t khit_n);
|
||||
|
||||
#define kv_pushp_ol(type, v, p) do { \
|
||||
if ((v).length == (v).size) { \
|
||||
(v).list = (type*)realloc((v).list, sizeof(type)*((v).size?((v).size<<1):(2))); \
|
||||
|
||||
+22
-11
@@ -712,19 +712,30 @@ inline uint32_t pop_trace_bp_f(asg16_v *res, uint32_t i, uint16_t *c, uint16_t *
|
||||
|
||||
inline int64_t pop_trace_bp_rev_f(asg16_v *res, int64_t i, uint16_t *c, uint16_t *bq, uint16_t *bt, uint32_t *len)
|
||||
{
|
||||
(*c) = (res->a[i]>>14); (*bq) = (*bt) = (uint16_t)-1;
|
||||
if((*c) == 2 || (*c) == 3) {
|
||||
(*bt) = ((res->a[i]>>12)&3);
|
||||
(*len) = (res->a[i]&(0xfff));
|
||||
} else if((*c) == 1) {
|
||||
(*bt) = ((res->a[i]>>12)&3);
|
||||
(*bq) = ((res->a[i]>>10)&3);
|
||||
(*len) = (res->a[i]&(0x3ff));
|
||||
} else {
|
||||
(*len) = (res->a[i]&(0x3fff));
|
||||
uint32_t sl = 1;
|
||||
(*c) = (*bq) = (*bt) = (uint16_t)-1; (*len) = 0;
|
||||
if(i >= ((int64_t)res->n)) {
|
||||
i = ((int64_t)res->n) - 1; sl = 0;
|
||||
}
|
||||
if(i >= 0) {
|
||||
(*c) = (res->a[i]>>14);
|
||||
(*bq) = (*bt) = (uint16_t)-1;
|
||||
if((*c) == 2 || (*c) == 3) {
|
||||
(*bt) = ((res->a[i]>>12)&3);
|
||||
(*len) = (res->a[i]&(0xfff));
|
||||
} else if((*c) == 1) {
|
||||
(*bt) = ((res->a[i]>>12)&3);
|
||||
(*bq) = ((res->a[i]>>10)&3);
|
||||
(*len) = (res->a[i]&(0x3ff));
|
||||
} else {
|
||||
(*len) = (res->a[i]&(0x3fff));
|
||||
}
|
||||
}
|
||||
if(sl == 0) {
|
||||
(*len) = 0; i = res->n;
|
||||
}
|
||||
|
||||
uint32_t sl; uint16_t sbq, sbt;
|
||||
uint16_t sbq, sbt;
|
||||
for (i--; (i >= 0) && ((*c) == (res->a[i]>>14)); i--) {
|
||||
sbq = sbt = (uint16_t)-1;
|
||||
if((*c) == 2 || (*c) == 3) {
|
||||
|
||||
+149
-18
@@ -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;
|
||||
@@ -23747,7 +23771,7 @@ void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sou
|
||||
sprintf(gfa_name, "%s.ovlp.reverse", output_file_name);
|
||||
write_ma_hit_ts(reverse_sources, RNF->total_reads, gfa_name);
|
||||
|
||||
if(asm_opt.dbg_bam) {
|
||||
/**if(asm_opt.dbg_bam)**/ {
|
||||
sprintf(gfa_name, "%s.rec", output_file_name);
|
||||
write_cc_v(&scb, gfa_name);
|
||||
}
|
||||
@@ -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);
|
||||
@@ -23794,7 +23818,7 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s
|
||||
return 0;
|
||||
}
|
||||
|
||||
if(asm_opt.dbg_bam) {
|
||||
/**if(asm_opt.dbg_bam)**/ {
|
||||
sprintf(gfa_name, "%s.rec", output_file_name);
|
||||
load_cc_v(&scb, gfa_name); ///write_ec_reads("lec.raw.fa", &scb, 0);
|
||||
}
|
||||
@@ -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);
|
||||
@@ -32495,6 +32559,11 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u
|
||||
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))
|
||||
{
|
||||
return 1;
|
||||
@@ -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);
|
||||
// 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,21 +40361,40 @@ 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,
|
||||
output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, cmk, 0);
|
||||
|
||||
+24
-7
@@ -125,6 +125,11 @@ void write_All_reads(All_reads* r, char* read_file_name)
|
||||
}
|
||||
}
|
||||
|
||||
mm = 3;
|
||||
fwrite(&mm, sizeof(mm), 1, fp);
|
||||
fwrite(&(r->tr[0]), sizeof(r->tr[0]), 1, fp);
|
||||
fwrite(&(r->tr[1]), sizeof(r->tr[1]), 1, fp);
|
||||
|
||||
free(index_name);
|
||||
fflush(fp);
|
||||
fclose(fp);
|
||||
@@ -235,11 +240,18 @@ int load_All_reads(All_reads* r, char* read_file_name)
|
||||
MALLOC(r->rsc[i], (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0));
|
||||
f_flag += fread(r->rsc[i], sizeof(uint8_t), (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0), fp);
|
||||
}
|
||||
|
||||
mm = 0;
|
||||
if(!feof(fp)) fread(&mm, sizeof(mm), 1, fp);
|
||||
}
|
||||
|
||||
if(mm == 3) {
|
||||
fread(&(r->tr[0]), sizeof(r->tr[0]), 1, fp);
|
||||
fread(&(r->tr[1]), sizeof(r->tr[1]), 1, fp);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
free(index_name);
|
||||
fclose(fp);
|
||||
fprintf(stderr, "Reads has been loaded.\n");
|
||||
@@ -283,12 +295,17 @@ uint8_t load_cc_v(cc_v* r, char* read_file_name)
|
||||
///typedef struct {size_t n, m; asg16_v *a; uint8_t *f; uint16_t *er; uint64_t bid;} cc_v;
|
||||
asg16_v *z; int f_flag = 0;
|
||||
uint64_t k, rn; uint32_t zn; f_flag += fread(&rn, sizeof(rn), 1, fp);
|
||||
r->n = r->m = rn; MALLOC(r->a, r->n);
|
||||
for (k = 0; k < r->n; k++) {
|
||||
z = &(r->a[k]);
|
||||
f_flag += fread(&zn, sizeof(zn), 1, fp);
|
||||
z->n = z->m = zn; MALLOC(z->a, z->n);
|
||||
f_flag += fread(z->a, sizeof((*(z->a))), zn, fp);
|
||||
r->n = r->m = rn;
|
||||
if(rn == 0) {
|
||||
free(r->a); r->a = NULL;
|
||||
} else {
|
||||
MALLOC(r->a, r->n);
|
||||
for (k = 0; k < r->n; k++) {
|
||||
z = &(r->a[k]);
|
||||
f_flag += fread(&zn, sizeof(zn), 1, fp);
|
||||
z->n = z->m = zn; MALLOC(z->a, z->n);
|
||||
f_flag += fread(z->a, sizeof((*(z->a))), zn, fp);
|
||||
}
|
||||
}
|
||||
|
||||
fflush(fp); fclose(fp);
|
||||
|
||||
+150
@@ -4014,6 +4014,156 @@ void get_pi_ec_chain(ha_abuf_t *ab, uint64_t rid, uint64_t rl, uint32_t tid, cha
|
||||
}
|
||||
|
||||
|
||||
uint64_t gen_srt_chain(ha_abuf_t *ab, ha_mz1_t *rfa, uint64_t *rfi, uint64_t rfn, ha_mz1_t *qra, uint64_t *qri, uint64_t qrn, uint32_t *high_occ, uint32_t *low_occ, Candidates_list *cl, uint32_t qrid)
|
||||
{
|
||||
uint64_t i, k, l, max_cnt = UINT32_MAX, min_cnt = 0, ri, qi, sn, x, rz, qz, n_a0, l0, k0, i0; anchor1_t *an; k_mer_hit *p;
|
||||
if(high_occ) {
|
||||
max_cnt = (*high_occ);
|
||||
if(max_cnt < 2) max_cnt = 2;
|
||||
}
|
||||
if(low_occ) {
|
||||
min_cnt = (*low_occ);
|
||||
if(min_cnt < 2) min_cnt = 2;
|
||||
}
|
||||
|
||||
radix_sort_anc64(rfi, rfi + rfn);
|
||||
radix_sort_anc64(qri, qri + qrn);
|
||||
i = ab->n_a = 0; n_a0 = l0 = i0 = 0; k0 = 1;
|
||||
for (k = 1, l = 0; k <= rfn; ++k) {
|
||||
if (k == rfn || (rfi[k]>>32) != (rfi[l]>>32)) {
|
||||
x = rfi[l]>>32;
|
||||
for (; i < qrn && (qri[i]>>32) < x; i++);
|
||||
if(ab->n_a < ab->m_a) {
|
||||
n_a0 = ab->n_a; l0 = l; i0 = i; k0 = k;
|
||||
}
|
||||
if(i < qrn && (qri[i]>>32) == x) {
|
||||
// sn = 0;
|
||||
// if(ab->n_a < ab->m_a) sn = ab->seed[l].n;///ha_pt_cnt(ha_idx, ra[l].x);
|
||||
for (qi = i; qi < qrn && (qri[qi]>>32) == x; qi++) {
|
||||
qz = (uint32_t)qri[qi];
|
||||
for (ri = l; ri < k; ri++) {
|
||||
rz = (uint32_t)rfi[ri];
|
||||
if(qra[qz].rev != rfa[rz].rev) continue;
|
||||
if(ab->n_a < ab->m_a) {
|
||||
an = &(ab->a[ab->n_a++]);
|
||||
an->other_off = qra[qz].pos;
|
||||
an->self_off = rfa[rz].pos;
|
||||
///an->cnt: cnt<<8|span
|
||||
an->cnt = ab->seed[rz].n; if(an->cnt > ((uint32_t)(0xffffffu))) an->cnt = 0xffffffu;
|
||||
an->cnt <<= 8; an->cnt |= ((rfa[rz].span <= ((uint32_t)(0xffu)))?rfa[rz].span:((uint32_t)(0xffu)));
|
||||
an->srt = (((uint64_t)(an->self_off))<<32)|((uint64_t)(an->other_off));
|
||||
} else {
|
||||
ab->n_a++;
|
||||
}
|
||||
}
|
||||
}
|
||||
i = qi;
|
||||
}
|
||||
l = k;
|
||||
}
|
||||
}
|
||||
|
||||
if (ab->n_a > ab->m_a) {
|
||||
ab->m_a = ab->n_a; REALLOC(ab->a, ab->m_a);
|
||||
k = k0; i = i0; l = l0; ab->n_a = n_a0;
|
||||
|
||||
for (; k <= rfn; ++k) {
|
||||
if (k == rfn || (rfi[k]>>32) != (rfi[l]>>32)) {
|
||||
x = rfi[l]>>32;
|
||||
for (; i < qrn && (qri[i]>>32) < x; i++);
|
||||
if(ab->n_a < ab->m_a) {
|
||||
n_a0 = ab->n_a; l0 = l; i0 = i; k0 = k;
|
||||
}
|
||||
if(i < qrn && (qri[i]>>32) == x) {
|
||||
// sn = 0;
|
||||
// if(ab->n_a < ab->m_a) sn = ab->seed[l].n;///ha_pt_cnt(ha_idx, ra[l].x);
|
||||
for (qi = i; qi < qrn && (qri[qi]>>32) == x; qi++) {
|
||||
qz = (uint32_t)qri[qi];
|
||||
for (ri = l; ri < k; ri++) {
|
||||
rz = (uint32_t)rfi[ri];
|
||||
if(qra[qz].rev != rfa[rz].rev) continue;
|
||||
if(ab->n_a < ab->m_a) {
|
||||
an = &(ab->a[ab->n_a++]);
|
||||
an->other_off = qra[qz].pos;
|
||||
an->self_off = rfa[rz].pos;
|
||||
///an->cnt: cnt<<8|span
|
||||
an->cnt = ab->seed[rz].n; if(an->cnt > ((uint32_t)(0xffffffu))) an->cnt = 0xffffffu;
|
||||
an->cnt <<= 8; an->cnt |= ((rfa[rz].span <= ((uint32_t)(0xffu)))?rfa[rz].span:((uint32_t)(0xffu)));
|
||||
an->srt = (((uint64_t)(an->self_off))<<32)|((uint64_t)(an->other_off));
|
||||
} else {
|
||||
ab->n_a++;
|
||||
}
|
||||
}
|
||||
}
|
||||
i = qi;
|
||||
}
|
||||
l = k;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// copy over to _cl_
|
||||
sn = ab->n_a + cl->length;
|
||||
if (sn > (uint64_t)cl->size) {
|
||||
cl->size = sn;
|
||||
REALLOC(cl->list, cl->size);
|
||||
}
|
||||
|
||||
radix_sort_ha_an1(ab->a, ab->a + ab->n_a);
|
||||
for (i = 0; i < ab->n_a; i++) {
|
||||
p = &cl->list[cl->length++];
|
||||
p->readID = qrid;
|
||||
p->strand = 0;
|
||||
p->offset = ab->a[i].other_off;
|
||||
p->self_offset = ab->a[i].self_off;
|
||||
|
||||
if(((ab->a[i].cnt>>8) < max_cnt) && ((ab->a[i].cnt>>8) > min_cnt)){
|
||||
p->cnt = 1;
|
||||
} else if((ab->a[i].cnt>>8) <= min_cnt) {
|
||||
p->cnt = 2;
|
||||
} else{
|
||||
p->cnt = 1 + (((ab->a[i].cnt>>8) + (max_cnt<<1) - 1)/(max_cnt<<1));
|
||||
p->cnt = pow(p->cnt, 1.1);
|
||||
}
|
||||
if(p->cnt > ((uint32_t)(0xffffffu))) p->cnt = 0xffffffu;
|
||||
p->cnt <<= 8; p->cnt |= (((uint32_t)(0xffu))&(ab->a[i].cnt));
|
||||
}
|
||||
// cl->length = ab->n_a;
|
||||
return ab->n_a;
|
||||
}
|
||||
|
||||
void gen_self_global_chain(ha_abuf_t *ab, Candidates_list *cl, uint32_t rid, uint64_t rl, uint32_t tid, char *ts, uint64_t tl, uint64_t mz_w, uint64_t mz_k, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, asg64_v *ix,
|
||||
uint8_t is_accurate, double bw_thres, int apend_be, uint32_t gen_off, int64_t enable_mcopy, double mcopy_rate, uint32_t mcopy_khit_cut, overlap_region_alloc *ores)
|
||||
{
|
||||
extern void *ha_flt_tab;
|
||||
uint64_t k, rn = ab->mz.n, m = cl->length;
|
||||
|
||||
mz1_ha_sketch(ts, tl, mz_w, mz_k, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL);
|
||||
|
||||
kv_resize(uint64_t, *ix, ab->mz.n);
|
||||
for (k = 0; k < rn; k++) {
|
||||
ix->a[k] = ab->mz.a[k].x;
|
||||
ix->a[k] <<= 32; ix->a[k] |= k;
|
||||
}
|
||||
|
||||
for (; k < ab->mz.n; k++) {
|
||||
ix->a[k] = ab->mz.a[k].x;
|
||||
ix->a[k] <<= 32; ix->a[k] |= (k-rn);
|
||||
}
|
||||
ix->n = ab->mz.n;
|
||||
|
||||
if(gen_srt_chain(ab, ab->mz.a, ix->a, rn, ab->mz.a + rn, ix->a + rn, ab->mz.n - rn, high_occ, low_occ, cl, tid)) {
|
||||
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);
|
||||
|
||||
m += lchain_qdp_mcopy_fast(cl, m, cl->length - m, m, &(cl->chainDP), ores, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres,
|
||||
rid, rl, tl, quick_check, apend_be, gen_off, enable_mcopy, mcopy_rate, mcopy_khit_cut, 1);
|
||||
cl->length = m;
|
||||
}
|
||||
ab->mz.n = rn;
|
||||
}
|
||||
|
||||
|
||||
int64_t ug_map_lchain(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, double bw_thres_sec,
|
||||
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, kv_u_trans_t *kov)
|
||||
|
||||
+1051
-89
File diff suppressed because it is too large
Load Diff
+178
-3
@@ -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;
|
||||
@@ -3053,44 +3200,57 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
|
||||
// fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n",
|
||||
// (int)Get_NAME_LENGTH(R_INF, 10819), Get_NAME(R_INF, 10819), 10819, is_contain_r((*rI), 10819));
|
||||
// debug_info_of_specfic_node("m64011_190830_220126/47516220/ccs", sg, rI, "beg-0");
|
||||
// debug_info_of_specfic_node("bcb40bcc-d9cf-48e6-88ee-47ac3dde22ff", sg, rI, "beg-0");
|
||||
// debug_info_of_specfic_node("c7ecbd6b-e09d-4042-93ac-2400839feaf6", sg, rI, "beg-0");
|
||||
|
||||
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);
|
||||
@@ -3106,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--");
|
||||
|
||||
/**
|
||||
@@ -3132,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;
|
||||
@@ -3146,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---");
|
||||
|
||||
@@ -3154,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);
|
||||
@@ -3183,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---");
|
||||
/**
|
||||
@@ -3197,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);
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -9,6 +9,7 @@
|
||||
#include "ksort.h"
|
||||
#include "htab.h"
|
||||
#include "Process_Read.h"
|
||||
#include "ecovlp.h"
|
||||
|
||||
#define YAK_COUNTER_BITS 12
|
||||
#define YAK_N_COUNTS (1<<YAK_COUNTER_BITS)
|
||||
@@ -1419,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);
|
||||
@@ -1478,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__);
|
||||
@@ -1486,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);
|
||||
@@ -1603,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);
|
||||
|
||||
@@ -1618,21 +1653,20 @@ 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, 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;
|
||||
|
||||
|
||||
if(is_load) {
|
||||
sprintf(gfa_name, "%s.r%lu.ht.bin", file_name, rr);
|
||||
fp = fopen(gfa_name, "r");
|
||||
@@ -1648,18 +1682,32 @@ uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
|
||||
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;
|
||||
free(gfa_name);
|
||||
return 0;
|
||||
}
|
||||
|
||||
sprintf(gfa_name, "%s.r%lu.rec", file_name, rr);
|
||||
if(!load_cc_v(rcc, gfa_name)) {
|
||||
destory_All_reads(r);
|
||||
ha_pt_destroy(*r_ha_idx); (*r_ha_idx) = NULL;
|
||||
ha_ft_destroy(*r_flt_tab); (*r_flt_tab) = NULL;
|
||||
destroy_cc_v(rcc);
|
||||
free(gfa_name);
|
||||
return 0;
|
||||
}
|
||||
|
||||
fprintf(stderr, "[M::%s] restarting from round %lu\n", __func__, rr0);
|
||||
} else {
|
||||
sprintf(gfa_name, "%s.r%lu.rec", file_name, rr);
|
||||
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");
|
||||
|
||||
@@ -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, 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);
|
||||
|
||||
@@ -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();
|
||||
|
||||
Reference in New Issue
Block a user