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

r852 -> hybrid correction

Browse Source
This commit is contained in:
chhylp123
2026-02-11 16:38:27 -05:00
parent ec9a8b222d
commit d6ba102d42
73 changed files with 190525 additions and 180074 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+41 -14
View File
@@ -883,16 +883,14 @@ static void worker_ec_save(void *data, long i, int tid)
void Output_corrected_reads()
{
long long i;
UC_Read g_read;
uint64_t i; UC_Read g_read;
init_UC_Read(&g_read);
char* gfa_name = (char*)malloc(strlen(asm_opt.output_file_name)+35);
sprintf(gfa_name, "%s.ec.fa", asm_opt.output_file_name);
FILE *output_file = fopen(gfa_name, "w");
free(gfa_name);
for (i = 0; i < (long long)R_INF.total_reads; i++)
{
for (i = 0; i < R_INF.total_reads; i++) {
recover_UC_Read(&g_read, &R_INF, i);
fwrite(">", 1, 1, output_file);
fwrite(Get_NAME(R_INF, i), 1, Get_NAME_LENGTH(R_INF, i), output_file);
@@ -906,7 +904,7 @@ void Output_corrected_reads()
void Output_corrected_fastq()
{
long long i; uint64_t k;
uint64_t i, k;
UC_Read g_read; asg8_v dv;
init_UC_Read(&g_read); kv_init(dv);
char* gfa_name = (char*)malloc(strlen(asm_opt.output_file_name)+35);
@@ -914,7 +912,7 @@ void Output_corrected_fastq()
FILE* fp = fopen(gfa_name, "w");
free(gfa_name);
for (i = 0; i < (long long)R_INF.total_reads; i++) {
for (i = 0; i < R_INF.tqn; i++) {
recover_UC_Read(&g_read, &R_INF, i);
fprintf(fp, "@%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
fprintf(fp, "%.*s\n", (int32_t)g_read.length, g_read.seq);
@@ -923,6 +921,17 @@ void Output_corrected_fastq()
for (k = 0; k < dv.n; k++) fprintf(fp, "%c", (char)(sc_tb[dv.a[k]] + 33 - 1));
fprintf(fp, "\n");
}
for (; i < R_INF.total_reads; i++) {
recover_UC_Read(&g_read, &R_INF, i);
fprintf(fp, "@%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
fprintf(fp, "%.*s\n", (int32_t)g_read.length, g_read.seq);
fprintf(fp, "+\n");
// retrive_bqual(&dv, NULL, i, -1, -1, 0, sc_bn);
// for (k = 0; k < dv.n; k++) fprintf(fp, "%c", (char)(sc_tb[dv.a[k]] + 33 - 1));
for (k = 0; k < (uint64_t)g_read.length; k++) fprintf(fp, "%c", (char)(3 + 33 - 1));
fprintf(fp, "\n");
}
destory_UC_Read(&g_read); kv_destroy(dv);
fclose(fp);
}
@@ -993,7 +1002,7 @@ void prt_dbg_rs(FILE *fp, Debug_reads* x, uint64_t round)
destory_UC_Read(&g_read);
}
void ha_ec(int64_t round, int num_pround, int des_idx, uint64_t *tot_b, uint64_t *tot_e)
void ha_ec(int64_t round, int num_pround, int des_idx, uint64_t *tot_b, uint64_t *tot_e, uint64_t w_tmp)
{
int hom_cov, het_cov, r_out = 0;
ha_flt_tab_hp = ha_idx_hp = NULL; (*tot_b) = (*tot_e) = 0;
@@ -1013,7 +1022,11 @@ void ha_ec(int64_t round, int num_pround, int des_idx, uint64_t *tot_b, uint64_t
het_cnt = NULL;
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);
if (r_out) {
write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
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);
// Output_corrected_fastq();
@@ -1955,8 +1968,13 @@ void ha_ec_ff(int renew_idx)
cal_ov_r(asm_opt.thread_num, R_INF.total_reads, renew_idx);
if(asm_opt.write_pos_idx) {
refresh_pt_idx(&ha_flt_tab, &ha_idx, NULL, &asm_opt, asm_opt.output_file_name, 1);
// write_pt_index(ha_flt_tab, ha_idx, NULL, &asm_opt, asm_opt.output_file_name);
} else {
ha_pt_destroy(ha_idx); ha_idx = NULL;
}
}
static void worker_ov_utg(void *data, long i, int tid)
{
@@ -2058,7 +2076,7 @@ int ha_assemble(void)
// debug_mc_gg_t(MC_NAME, 0, 0);
// quick_debug_phasing(MC_NAME);
extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
int r, hom_cov = -1, ovlp_loaded = 0; uint64_t tot_b, tot_e;
int r, r0 = -1, hom_cov = -1, ovlp_loaded = 0; uint64_t tot_b, tot_e;
if (asm_opt.load_index_from_disk && 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());
@@ -2076,20 +2094,29 @@ 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);
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)) {
load_ct_index(&ha_ct_table, asm_opt.output_file_name); r0 = r;
break;
}
}
if(r < 0) r = 0;
}
// construct hash table for high occurrence k-mers
if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL)
{
if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL) {
ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0, 0);
ha_opt_update_cov(&asm_opt, hom_cov);
}
// error correction
assert(asm_opt.number_of_round > 0);
for (r = ha_idx?asm_opt.number_of_round-1:0; r < asm_opt.number_of_round; ++r) {
for (; r < asm_opt.number_of_round; ++r) {
ha_opt_reset_to_round(&asm_opt, r); // this update asm_opt.roundID and a few other fields
tot_b = tot_e = 0;
// ha_overlap_and_correct(r);
ha_ec(r, asm_opt.number_of_pround, (r<asm_opt.number_of_round-1)?1:0, &tot_b, &tot_e);
ha_ec(r, asm_opt.number_of_pround, (r<asm_opt.number_of_round-1)?1:0, &tot_b, &tot_e, ((r > r0) && (asm_opt.restart))?1:0);
fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(),
yak_cpu_usage(), yak_peakrss_in_gb(), r + 1);
fprintf(stderr, "[M::%s] # bases: %lu; # corrected bases: %lu\n", __func__, tot_b, tot_e);
@@ -2108,7 +2135,7 @@ int ha_assemble(void)
fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb());
// 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);
ha_ft_destroy(ha_flt_tab);
if(!(asm_opt.write_pos_idx)) ha_ft_destroy(ha_flt_tab);
if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
ha_triobin(&asm_opt);
CommandLines.cpp
+85 -12
View File
@@ -81,6 +81,13 @@ static ko_longopt_t long_options[] = {
{ "ul-m", ko_required_argument, 363},
{ "rl-cut", ko_required_argument, 364},
{ "sc-cut", ko_required_argument, 365},
{ "hf", ko_required_argument, 366},
{ "cb", ko_required_argument, 367},
{ "gpath", ko_no_argument, 368},
{ "het-cov", ko_required_argument, 369},
{ "resume", ko_no_argument, 370},
{ "flt-kocc", ko_required_argument, 371},
{ "chn-occ", ko_required_argument, 372},
// { "path-round", ko_required_argument, 348},
{ 0, 0, 0 }
};
@@ -107,7 +114,9 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " -k INT k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
fprintf(stderr, " -w INT minimizer window size [%d]\n", asm_opt->mz_win);
fprintf(stderr, " -f INT number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift);
fprintf(stderr, " -D FLOAT drop k-mers occurring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
fprintf(stderr, " -D FLOAT drop k-mers occurring >FLOAT*coverage times [%.1f]; work with --flt-kocc or -N\n", asm_opt->high_factor);
fprintf(stderr, " --flt-kocc INT\n");
fprintf(stderr, " drop k-mers occurring >max(-D*coverage,--flt-kocc) times [%ld]\n", asm_opt->hf_cutoff);
fprintf(stderr, " -N INT consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
fprintf(stderr, " -r INT round of correction [%d]\n", asm_opt->number_of_round);
fprintf(stderr, " -z INT length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
@@ -115,6 +124,13 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " employ k-mers occurring <INT times to rescue repetitive overlaps [%d]\n", asm_opt->max_kmer_cnt);
fprintf(stderr, " --hg-size INT(k, m or g)\n");
fprintf(stderr, " estimated haploid genome size used for inferring read coverage [auto]\n");
fprintf(stderr, " --resume resume from the previously incomplete assembly [%ld]\n", asm_opt->restart);
fprintf(stderr, " --het-cov INT\n");
fprintf(stderr, " heterozygous read coverage [auto]; used for error correction and assembly; manual value overrides auto\n");
fprintf(stderr, " --hom-cov INT\n");
fprintf(stderr, " homozygous read coverage [auto]; used for error correction and assembly; manual value overrides auto\n");
fprintf(stderr, " --chn-occ INT\n");
fprintf(stderr, " discard overlaps supported by <INT minimizers [%ld]\n", asm_opt->chn_occ);
fprintf(stderr, " Assembly:\n");
fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round);
fprintf(stderr, " -m INT pop bubbles of <INT in size in contig graphs [%lld]\n", asm_opt->large_pop_bubble_size);
@@ -126,8 +142,6 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " -u post-join step for contigs which may improve N50; 0 to disable; 1 to enable\n");
fprintf(stderr, " [%u] and [%u] in default for the UL+HiFi assembly and the HiFi assembly, respectively\n",
asm_opt->ul_pst_join, asm_opt->hifi_pst_join);
fprintf(stderr, " --hom-cov INT\n");
fprintf(stderr, " homozygous read coverage [auto]\n");
fprintf(stderr, " --lowQ INT\n");
fprintf(stderr, " output contig regions with >=INT%% inconsistency in BED format; 0 to disable [%d]\n", asm_opt->bed_inconsist_rate);
fprintf(stderr, " --b-cov INT\n");
@@ -140,6 +154,7 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " --primary output a primary assembly and an alternate assembly\n");
fprintf(stderr, " --ctg-n INT\n");
fprintf(stderr, " remove tip contigs composed of <=INT reads [%d]\n", asm_opt->max_contig_tip);
fprintf(stderr, " --gpath output the corresponding path of each contig (p_ctg) within the assembly graph (d_utg.noseq.gfa)\n");
// fprintf(stderr, " --pri-range INT1[,INT2]\n");
@@ -223,17 +238,18 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " --telo-s INT\n");
fprintf(stderr, " min score for telomere reads [%ld]\n", asm_opt->telo_mic_sc);
fprintf(stderr, " ONT simplex assembly (beta):\n");
fprintf(stderr, " --ont assemble ONT simplex reads in fastq format\n");
fprintf(stderr, " ONT Simplex assembly (beta):\n");
fprintf(stderr, " --ont assemble ONT Simplex reads in fastq format\n");
// fprintf(stderr, " --sc-n consider base qual value for assembly\n");
fprintf(stderr, " --chem-c INT\n");
fprintf(stderr, " detect chimeric reads with <=INT other reads support [%lu]\n", asm_opt->chemical_cov);
fprintf(stderr, " --chem-f INT\n");
fprintf(stderr, " length of flanking regions for chimeric read detection [%lu]\n", asm_opt->chemical_flank);
fprintf(stderr, " --rl-cut INT\n");
fprintf(stderr, " filter out ONT simplex reads shorter than <INT> for assembly [%ld]\n", asm_opt->rl_cut);
fprintf(stderr, " filter out ONT Simplex reads shorter than <INT> for assembly [%ld]\n", asm_opt->rl_cut);
fprintf(stderr, " --sc-cut INT\n");
fprintf(stderr, " filter out ONT simplex reads with a mean base quality score below <INT> [%ld]\n", asm_opt->sc_cut);
fprintf(stderr, " filter out ONT Simplex reads with a mean base quality score below <INT> [%ld]\n", asm_opt->sc_cut);
fprintf(stderr, " --hf FILEs file names of HiFi reads\n");
fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
@@ -254,6 +270,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->hic_reads[0] = NULL;
asm_opt->hic_reads[1] = NULL;
asm_opt->fn_bin_poy = NULL;
asm_opt->fn_chr_bin = NULL;
asm_opt->ar = NULL;
asm_opt->thread_num = 1;
asm_opt->k_mer_length = 51;
@@ -270,6 +287,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->max_kmer_cnt = 2000;
asm_opt->high_factor = 5.0;
asm_opt->max_ov_diff_ec = 0.04;
asm_opt->max_ov_diff_ec_sec = 0.04;
asm_opt->max_ov_diff_final = 0.03;
asm_opt->hom_cov = 20;
asm_opt->het_cov = -1024;
@@ -373,6 +391,30 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->rl_cut = 1000;
asm_opt->sc_cut = 10;
asm_opt->hf = NULL;
asm_opt->gpath = 0;
asm_opt->hf_rate = 4;
asm_opt->ont_rate = 1;///must be 1 or 0
asm_opt->hf_rate_max = 4;
asm_opt->het_cov_set = -1;
asm_opt->restart = 0;
asm_opt->hf_cutoff = -1;
asm_opt->write_pos_idx = 1;
asm_opt->hom_cov_0 = -1;
asm_opt->het_cov_0 = -1;
asm_opt->max_n_chain_0 = -1;
asm_opt->hmo_cov_ss = -1;
asm_opt->het_cov_ss = -1;
asm_opt->chn_occ = 2;
}
void destory_enzyme(enzyme* f)
@@ -752,6 +794,12 @@ int check_option(hifiasm_opt_t* asm_opt)
}
}
if((asm_opt->hf) && (!(asm_opt->is_ont))) {
fprintf(stderr, "[ERROR] [--hf] must work with [--ont]\n");
return 0;
}
return 1;
}
@@ -920,12 +968,14 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 306) asm_opt->max_ov_diff_final = atof(opt.arg);
else if (c == 307) asm_opt->extract_list = opt.arg;
else if (c == 308) asm_opt->extract_iter = atoi(opt.arg);
else if (c == 309)
{
asm_opt->hom_global_coverage = atoi(opt.arg);
else if (c == 309) {
asm_opt->hom_global_coverage = asm_opt->hmo_cov_ss = atoi(opt.arg);
asm_opt->hom_global_coverage_set = 1;
if(asm_opt->hmo_cov_ss <= 0) {
fprintf(stderr, "[ERROR] homozygous read coverage should be > 0 (--hom-cov)");
return 1;
}
else if (c == 310)
} else if (c == 310)
{
char* s = NULL;
asm_opt->recover_atg_cov_min = strtol(opt.arg, &s, 10);
@@ -1003,7 +1053,30 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
asm_opt->rl_cut = atol(opt.arg);
} else if (c == 365) {
asm_opt->sc_cut = atol(opt.arg);
} else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap
} else if (c == 366) {
get_hic_enzymes(opt.arg, &(asm_opt->hf), 0);
} else if (c == 367) {
asm_opt->fn_chr_bin = opt.arg;
} else if (c == 368) {
asm_opt->gpath = 1;
} else if (c == 369) {
asm_opt->het_cov_set = asm_opt->het_cov_ss = atoi(opt.arg);
if(asm_opt->het_cov_ss <= 0) {
fprintf(stderr, "[ERROR] heterozygous read coverage should be > 0 (--het-cov)");
return 1;
}
} else if (c == 370) {
asm_opt->restart = 1;
} else if (c == 371) {
asm_opt->hf_cutoff = atoi(opt.arg);
} else if (c == 372) {
asm_opt->chn_occ = atoi(opt.arg);
if(asm_opt->chn_occ <= 0) {
fprintf(stderr, "[ERROR] chain cutoff should be > 0 (--chn-occ)");
return 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);
}
else if (c == 's') asm_opt->purge_simi_rate_l2 = asm_opt->purge_simi_rate_l3 = atof(opt.arg);
CommandLines.h
+23 -1
View File
@@ -5,7 +5,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.25.0-r726"
#define HA_VERSION "0.25.0-r852"
#define VERBOSE 0
@@ -41,10 +41,12 @@ typedef struct {
char *fn_bin_yak[2];
char *fn_bin_list[2];
char *fn_bin_poy;
char *fn_chr_bin;
char *extract_list;
enzyme *hic_reads[2];
enzyme *hic_enzymes;
enzyme *ar;
enzyme *hf;
enzyme *sec_in;
int extract_iter;
int thread_num;
@@ -63,6 +65,7 @@ typedef struct {
int max_kmer_cnt;
double high_factor; // coverage cutoff set to high_factor*hom_cov
double max_ov_diff_ec;
double max_ov_diff_ec_sec;
double max_ov_diff_final;
int hom_cov;
int het_cov;
@@ -169,7 +172,26 @@ typedef struct {
int64_t rl_cut;
int64_t sc_cut;
uint8_t gpath;
uint64_t hf_rate;///cannot be larger than 128?
uint64_t hf_rate_max;///cannot be larger than 128?
uint64_t ont_rate;///cannot be 0, should be 1 in anyway
int64_t het_cov_set;
int64_t restart;
int64_t hf_cutoff;
uint8_t write_pos_idx;
int hom_cov_0;
int het_cov_0;
int max_n_chain_0; // fall-back max number of chains to consider
int64_t hmo_cov_ss;
int64_t het_cov_ss;
int64_t chn_occ;
} hifiasm_opt_t;
extern hifiasm_opt_t asm_opt;
Correct.cpp
+8706 -360
View File
File diff suppressed because it is too large Load Diff
Correct.h
+69 -3
View File
@@ -1387,18 +1387,71 @@ typedef struct {
int64_t k, q[2], t[2], cq[2], ct[2], ci[2], werr, werr0, cerr;
int64_t qoff, f, toff, coff, cur_qoff;
} rtrace_iter;
typedef struct
{
overlap_region_alloc *ol;
Candidates_list *cl;
All_reads *rref;
UC_Read *qu;
UC_Read *tu;
bit_extz_t *exz;
overlap_region *aux_o, *rse_o;
double e_rate[2];
int64_t wl[2];
int64_t rid;
int64_t khit;
int64_t move_gap;
asg16_v *buf;
asg64_v *srt;
asg8_v *hpz;
ma_hit_t_alloc *in;
int8_t chem_drop[2];
double align_gap_rate[2];
int64_t align_gap_max[2];
uint64_t sec_aln_win;
uint64_t sec_aln_cov;
double sec_aln_err_rate;
double sec_aln_max;
asg64_v *kp;
asg32_v *v32;
asg64_v *bp;
ha_abuf_t *ab;
uint64_t max_n_chain;
uint64_t max_n_chain_f;
uint64_t chain_cutoff;
uint64_t ave_cov_min;
uint64_t ocw;
uint64_t t_cut;
} gen_hc_aln_t;
int64_t get_rid_backward_cigar_err(rtrace_iter *it, ul_ov_t *aln, kv_rtrace_t *trace, rtrace_t *tc,
const ul_idx_t *uref, char* qstr, UC_Read *tu, overlap_region_alloc *ol, overlap_region *o,
bit_extz_t *exz, double e_rate, int64_t qs);
void gen_hc_r_alin(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);
void 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);
void gen_hc_r_alin_adp_smp(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 *sp, uint64_t ocw, uint8_t *hpf, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match);
void gen_hc_r_alin_adv_adp_smp(gen_hc_aln_t *ez, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, uint8_t set_match);
void pp_chn_a(overlap_region *z, Candidates_list *cl, uint8_t is_raw);
void gen_hc_r_alin(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_flt(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, overlap_region *aux_b, 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, asg64_v *sp, uint64_t ocw, uint8_t *hpf, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min);
void gen_hc_r_alin_adp(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, overlap_region *aux_b, 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, asg64_v *sp, uint64_t ocw, uint8_t *hpf, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min);
void gen_hc_r_alin_adv(gen_hc_aln_t *ez);
void 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);
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);
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);
void get_wqual(uint64_t zid, uint64_t zpos, uint64_t zrev, asg8_v *v, uint8_t *va, uint64_t scw, uint64_t *tqual, uint64_t *wqual);
void gen_reseed_re(overlap_region_alloc *ol, Candidates_list *cl, overlap_region *aux_o, overlap_region *rse_o, All_reads *rref, UC_Read* qu, UC_Read *tu, bit_extz_t *exz, kv_ul_ov_t *c_idx, asg64_v *idx, asg64_v *res, int64_t bd, int64_t mzw, int64_t kl, int64_t rid, double err_h, double err_l, asg16_v *b16, uint64_t tqn, uint8_t *hpf);
#define ovlp_id(x) ((x).tn)
@@ -1410,11 +1463,24 @@ void get_wqual(uint64_t zid, uint64_t zpos, uint64_t zrev, asg8_v *v, uint8_t *v
#define ovlp_cur_ylen(x) ((x).te)
#define ovlp_cur_coff(x) ((x).qe)
#define ovlp_bd(x) ((x).sec)
#define ovlp_hf(x) ((x).el)
#define HPC_PL 12
#define HPC_RR 4
#define HPC_CC 2
#define HPC_RR_Q 5
#define HPC_CC_Q 3
#define HC_MF_R 0.5
#define HC_AV_MIN 0.7
// #define FORCE_CUT 1
// 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, asg64_v *srt, ma_hit_t_alloc *in, 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
// uint64_t gen_hc_r_alin_ea_hybrid(overlap_region_alloc* ol, uint64_t bi, uint64_t bn, uint64_t tk, 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 ec_filter, ma_hit_t_alloc *in, 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)
#endif
Hash_Table.cpp
+58 -8
View File
@@ -1749,8 +1749,7 @@ uint64_t lchain_qdp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, o
return cL;
}
void push_ovlp_chain_qgen(overlap_region* o, uint32_t xid, int64_t xl, int64_t yl, int64_t sc,
k_mer_hit *beg, k_mer_hit *end)
void push_ovlp_chain_qgen(overlap_region* o, uint32_t xid, int64_t xl, int64_t yl, int64_t sc, k_mer_hit *beg, k_mer_hit *end)
{
int64_t xr, yr;
o->x_id = xid; o->y_id = beg->readID;
@@ -2108,9 +2107,9 @@ uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n,
t = dp->tmp; f = dp->score; p = dp->pre; ii = dp->occ;
a = cl->list + a_idx; des = cl->list + des_idx;
// if(a_n && a[0].readID == 0) {
// fprintf(stderr, "---[M::%s::utg%.6dl::%c]\n",
// __func__, (int32_t)a[0].readID+1, "+-"[a[0].strand]);
// if(a_n && (a[0].readID == 27105 || a[0].readID == 7603)) {///r833
// fprintf(stderr, "---[M::%s::rid->%u::%c]\ta_n::%ld\n",
// __func__, a[0].readID, "+-"[a[0].strand], a_n);
// }
if(quick_check) {
quick_ck_lchain(a, a_n, xl, yl, chn_pen_gap, chn_pen_skip, bw_rate, p, t, f, ii, &plus, &msc, &msc_i, &movl, &si, &ei);
@@ -2168,9 +2167,9 @@ uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n,
}
if(f[i] < plus) plus = f[i];
ii[i] = 0;///for mcopy, not here
// if(a_n && a[0].readID == 4412344) {
// fprintf(stderr, "i::%ld[M::%s::%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].strand],
// 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);
// }
}
@@ -2604,6 +2603,57 @@ uint64_t lchain_simple(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp
return cL;
}
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)
{
if(a_n <= 0) return 0;
int64_t *p, *t, max_f, n_skip, st, max_j, sc, msc, msc_i;
int32_t *f; int64_t i, j, cL = 0;
resize_Chain_Data(dp, a_n, NULL);
t = dp->tmp; f = dp->score; p = dp->pre; msc = msc_i = -1;
memset(t, 0, (a_n*sizeof((*t))));
f[0]=a[0].cnt; p[0]=-1; msc = f[0]; msc_i = 0;
for (i = 1, st = 0; i < a_n; ++i) {
max_f = INT32_MIN; n_skip = 0; max_j = -1;
if ((i-st) > max_iter) st = i-max_iter;
///[st, i-2]
for (j=i-1; j >= st; --j) {
if((a[i].self_offset > a[j].self_offset)&&(a[i].offset > a[j].offset)) {
sc = f[j]+a[i].cnt;
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;
}
}
f[i] = max_f; p[i] = max_j;
if(f[i] > msc) {
msc = f[i]; msc_i = i;
}
}
///a[] has been sorted by self_offset
i = msc_i;
cL = 0;
while (i >= 0) {
t[cL++] = 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(des) {
for (i = 0; i < cL; i++) des[i] = a[t[i]];
}
return cL;
}
inline int64_t hit_long_gap(k_mer_hit *a, k_mer_hit *b, int64_t max_lgap, double small_bw_rate, int64_t min_small_bw)
{
int64_t dq, dr, dd, dm;
Hash_Table.h
+5
View File
@@ -17,6 +17,8 @@
#define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25
#define THRESHOLD 15
#define OVERLAP_THRESHOLD_HIFI_FILTER 0.9
#define OVERLAP_THRESHOLD_HIFI_FF_FILTER 0.6
#define OVERLAP_THRESHOLD_HIFI_FF_DE_FILTER 0.5
#define OVERLAP_THRESHOLD_NOSI_FILTER 0.7
#define OVERLAP_THRESHOLD_FILTER_HPC 0.75
#define HIGH_HET_OVERLAP_THRESHOLD_FILTER 0.3
@@ -234,6 +236,9 @@ uint64_t lchain_qdp_fix(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_s
int64_t left_fix, int64_t right_fix);
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);
void push_ovlp_chain_qgen(overlap_region* o, uint32_t xid, int64_t xl, int64_t yl, int64_t sc, k_mer_hit *beg, k_mer_hit *end);
uint64_t lchain_qdp_mcopy(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,
Levenshtein_distance.h
+12 -2
View File
@@ -796,7 +796,7 @@ inline uint32_t cigar_check(char *pstr, char *tstr, bit_extz_t *ez)
if(c == 0) {
for (k=0;(k<cl)&&(pstr[pi]==tstr[ti]);k++,pi++,ti++);
if(k!=cl) {
fprintf(stderr, "ERROR-d-0, pi::%d, ti::%d\n", pi, ti);
fprintf(stderr, "ERROR-d-0, pi::%d, ti::%d, ci::%u\n", pi, ti, ci);
return 0;
}
} else {
@@ -804,7 +804,7 @@ inline uint32_t cigar_check(char *pstr, char *tstr, bit_extz_t *ez)
if(c == 1) {
for (k=0;(k<cl)&&(pstr[pi]!=tstr[ti]);k++,pi++,ti++);
if(k!=cl) {
fprintf(stderr, "ERROR-d-1, pi::%d, ti::%d\n", pi, ti);
fprintf(stderr, "ERROR-d-1, pi::%d, ti::%d, ci::%u\n", pi, ti, ci);
return 0;
}
} else if(c == 2) {///more p
@@ -814,10 +814,20 @@ inline uint32_t cigar_check(char *pstr, char *tstr, bit_extz_t *ez)
}
}
}
if(err != ez->err) {
fprintf(stderr, "ERROR-err\n");
return 0;
}
if(pi != ez->pe + 1) {
fprintf(stderr, "ERROR-pi\n");
return 0;
}
if(ti != ez->te + 1) {
fprintf(stderr, "ERROR-ti\n");
return 0;
}
return 1;
}
Overlaps.cpp
+554 -29
View File
@@ -75,7 +75,7 @@ KRADIX_SORT_INIT(ha_mzl_t_srt1, ha_mzl_t, ha_mzl_t_key, member_size(ha_mzl_t, x)
KSORT_INIT_GENERIC(uint32_t)
void reduce_hamming_error_adv(ma_ug_t *iug, asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub);
int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub, uint32_t max_ext);
void print_vw_edge(asg_t *sg, uint32_t vid, uint32_t wid, const char *cmd);
void output_trio_graph_joint(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long tipsLen, float tip_drop_ratio,
@@ -8966,7 +8966,7 @@ add_unitig:
if((ct != FATHER) && (ct != MOTHER)) ct = AMBIGU;
pt = ct;
fn0 = fn; mn0 = mn; uidx.n = fn = mn = 0;
if(ct == FATHER) fn++; if(ct == MOTHER) mn++;
if(ct == FATHER) {fn++;} if(ct == MOTHER) {mn++;}
for (k = 1, l = 0; k <= kdq_size(q); k++) {
st = 0; ct = AMBIGU;
if(k == kdq_size(q)) {
@@ -8990,7 +8990,7 @@ add_unitig:
fn = mn = 0; l = k;
}
if(ct != AMBIGU) pt = ct;
if(ct == FATHER) fn++; if(ct == MOTHER) mn++;
if(ct == FATHER) {fn++;} if(ct == MOTHER) {mn++;}
}
@@ -9025,7 +9025,7 @@ add_unitig:
p->a[i] = kdq_at(q, i);
ct = R_INF.trio_flag[p->a[i]>>33];
if((ct != FATHER) && (ct != MOTHER)) ct = AMBIGU;
if(ct == FATHER) fn1++; if(ct == MOTHER) mn1++;
if(ct == FATHER) {fn1++;} if(ct == MOTHER) {mn1++;}
}
} else {
p->s = 0; p->len = 0; p->circ = 0;
@@ -9038,13 +9038,13 @@ add_unitig:
p->a[z] = kdq_at(q, i); p->len += (uint32_t)p->a[z];
ct = R_INF.trio_flag[p->a[z]>>33];
if((ct != FATHER) && (ct != MOTHER)) ct = AMBIGU;
if(ct == FATHER) fn1++; if(ct == MOTHER) mn1++;
if(ct == FATHER) {fn1++;} if(ct == MOTHER) {mn1++;}
}
p->a[z] = kdq_at(q, i); p->a[z] >>= 32; p->a[z] <<= 32;
p->a[z] |= g->seq[p->a[z]>>33].len; p->len += (uint32_t)p->a[z];
ct = R_INF.trio_flag[p->a[z]>>33];
if((ct != FATHER) && (ct != MOTHER)) ct = AMBIGU;
if(ct == FATHER) fn1++; if(ct == MOTHER) mn1++;
if(ct == FATHER) {fn1++;} if(ct == MOTHER) {mn1++;}
}
fn = mn = 0; l = k;
if(k < uidx.n) {
@@ -10755,7 +10755,7 @@ uint32_t get_ug_coverage(ma_utg_t* u, asg_t* read_g, const ma_sub_t* coverage_cu
uint32_t k, j, rId, tn, is_Unitig;
long long R_bases = 0, C_bases = 0;
ma_hit_t *h;
if(rR) (*rR) = 0; if(rC) (*rC) = 0;
if(rR) {(*rR) = 0;} if(rC) {(*rC) = 0;}
if(u->m == 0) return 0;
for (k = 0; k < u->n; k++)
@@ -10795,7 +10795,7 @@ uint32_t get_ug_coverage(ma_utg_t* u, asg_t* read_g, const ma_sub_t* coverage_cu
r_flag[rId] = 0;
}
if(rR) (*rR) = R_bases; if(rC) (*rC) = C_bases;
if(rR) {(*rR) = R_bases;} if(rC) {(*rC) = C_bases;}
return R_bases == 0? 0:C_bases/R_bases;
}
@@ -11101,6 +11101,64 @@ void ma_scg_print(const kvect_sec_t *sc, All_reads *RNF, asg_t* read_g, const ma
free(primary_flag);
}
void ma_scg_gpath_print(const kvect_sec_t *sc, const char* prefix, scaf_res_t *gpt, ma_ug_t *gpt_u, FILE *fp)
{
uc_block_t *a; uint64_t i, j, k, zk, tl, scl, a_n; sec_t *scp; ul_vec_t *idx;
char name[32]; ma_utg_t *z = NULL;
for (i = 0; i < sc->n; ++i) { // the Segment lines in GFA
scp = &(sc->a[i]);
///debug
// prt_scaf_stats(scp, sc->ctg, prefix, i);
sprintf(name, "%s%.6lu%c", prefix, i + 1, "lc"[scp->is_c]);
for (j = scl = 0; j < scp->n; j++) {
z = &(sc->ctg->u.a[((uint32_t)scp->a[j])>>1]);
scl += z->len; scl += (scp->a[j]>>32);
}
for (j = tl = 0; j < scp->n; j++) {
z = &(sc->ctg->u.a[((uint32_t)scp->a[j])>>1]); ///each contig
idx = &(gpt->a[((uint32_t)scp->a[j])>>1]); ///alignment path
if(!(((uint32_t)scp->a[j])&1)) {///forward
for (k = 0; k < idx->bb.n; k++) {
a = idx->bb.a + idx->bb.a[k].ts; a_n = idx->bb.a[k].te - idx->bb.a[k].ts;
for (zk = 0; zk < a_n; zk++) {///alignment
fprintf(fp, "%s\t%lu\t%lu\t%lu\t%c\tutg%.6u%c\t%u\t%u\t%u\n", name, scl, tl + a[zk].qs, tl + a[zk].qe, "+-"[a[zk].rev],
(a[zk].hid)+1, "lc"[gpt_u->u.a[a[zk].hid].circ], gpt_u->u.a[a[zk].hid].len, a[zk].ts, a[zk].te);
}
}
} else {///revrse
for (k = 0; k < idx->bb.n; k++) {
a = idx->bb.a + idx->bb.a[idx->bb.n-k-1].ts; a_n = idx->bb.a[idx->bb.n-k-1].te - idx->bb.a[idx->bb.n-k-1].ts;
for (zk = 0; zk < a_n; zk++) {
fprintf(fp, "%s\t%lu\t%lu\t%lu\t%c\tutg%.6u%c\t%u\t%u\t%u\n", name, scl, tl + z->len - a[a_n-zk-1].qe, tl + z->len - a[a_n-zk-1].qs, "+-"[a[a_n-zk-1].rev^1],
(a[a_n-zk-1].hid)+1, "lc"[gpt_u->u.a[a[a_n-zk-1].hid].circ], gpt_u->u.a[a[a_n-zk-1].hid].len, a[a_n-zk-1].ts, a[a_n-zk-1].te);
}
}
}
tl += z->len; tl += (scp->a[j]>>32);
}
}
}
void ma_ug_gpath_print(const char* prefix, scaf_res_t *gpt, ma_ug_t *ctg_u, ma_ug_t *ref_u, FILE *fp)
{
uc_block_t *a; uint64_t i, k, zk, a_n; ul_vec_t *idx;
char name[32]; ma_utg_t *z = NULL;
for (i = 0; i < ctg_u->u.n; ++i) { // the Segment lines in GFA
///debug
// prt_scaf_stats(scp, sc->ctg, prefix, i);
z = &(ctg_u->u.a[i]); idx = &(gpt->a[i]); ///alignment path
sprintf(name, "%s%.6lu%c", prefix, i + 1, "lc"[z->circ]);
for (k = 0; k < idx->bb.n; k++) {
a = idx->bb.a + idx->bb.a[k].ts; a_n = idx->bb.a[k].te - idx->bb.a[k].ts;
for (zk = 0; zk < a_n; zk++) {///alignment
fprintf(fp, "%s\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\n", name, z->len, a[zk].qs, a[zk].qe, "+-"[a[zk].rev],
(a[zk].hid)+1, "lc"[ref_u->u.a[a[zk].hid].circ], ref_u->u.a[a[zk].hid].len, a[zk].ts, a[zk].te);
}
}
}
}
int asg_cut_internal(asg_t *g, int max_ext)
{
asg64_v a = {0,0,0};
@@ -14679,6 +14737,7 @@ void debug_hapS(uint32_t *hapS, uint32_t rn)
free(occ); free(freq);
}
}
void output_poly_trio(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, long long tipsLen, float tip_drop_ratio, long long stops_threshold,
R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int gap_fuzz, int is_bench,
@@ -14697,6 +14756,37 @@ bub_label_t* b_mask_t, uint32_t hapN)
free(fp); free(hapS);
}
void output_chr_bin_trio(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, long long tipsLen, float tip_drop_ratio, long long stops_threshold,
R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int gap_fuzz, int is_bench,
bub_label_t* b_mask_t)
{
uint32_t i, k, idx_n = 0, na; uint8_t *idx = NULL;
uint32_t *hapS = ha_charbin_list(&asm_opt, &idx, &idx_n);
fprintf(stderr, "[M::%s] # reads: %u\n", __func__, sg->n_seq);
// debug_hapS(hapS, sg->n_seq);
char *fp = NULL; MALLOC(fp, 100);
for (i = 0; i < idx_n; i++){
if(!idx[i]) continue;
for (k = na = 0; k < sg->n_seq; k++) {
if(R_INF.trio_flag[k] == DROP) continue;
R_INF.trio_flag[k] = AMBIGU;
if(hapS[k] == ((uint32_t)-1)) continue;
if(hapS[k] == i) {
R_INF.trio_flag[k] = FATHER; na++;
} else {
R_INF.trio_flag[k] = MOTHER;
}
}
if(!na) continue;
sprintf(fp, "hap%u", i);
fprintf(stderr, "[M::%s] # %s reads: %u\n", __func__, fp, na);
output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, reverse_sources, tipsLen, tip_drop_ratio,
stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, gap_fuzz, is_bench, b_mask_t, fp, NULL, NULL);
}
free(fp); free(hapS); free(idx);
}
uint32_t test_dbug(ma_ug_t* ug, FILE* fp)
{
uint32_t f_flag = 0, t, i, r_flag = 0;
@@ -15166,7 +15256,7 @@ void deep_clean_u_trans(kv_u_trans_t *ta, u_trans_t *mz, ma_ug_t *ug, asg64_v *s
mz->occ = cz->occ = 0; return;
} else {
a = u_trans_a(*ta, mz->qn); n = u_trans_n(*ta, mz->qn);
for (k = 0; (k < n) && (a[k].tn != mz->tn); k++); assert(k < n);
for (k = 0; (k < n) && (a[k].tn != mz->tn); k++){;} assert(k < n);
if(!test_arc_rm(ug, a, n, k, ug->u.a[mz->qn].len, srt, len_cut, rate_cut)) {
mz->occ = cz->occ = 0; return;
}
@@ -15562,7 +15652,7 @@ void clean_u_trans_t_idx_filter_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g
sl.cu = gen_u_trans_cluster(ta);
kt_for(sl.n_thread, worker_for_trans_sec_cut, &sl, sl.cu->idx.n);
free(sl.cu->idx.a); free(sl.cu->z.a); free(sl.cu);
for (k = 0; k < sl.n_thread; k++) free(sl.srt[k].a); free(sl.srt);
for (k = 0; k < sl.n_thread; k++) {free(sl.srt[k].a);} free(sl.srt);
for (k = st = 0; k < ta->n; k++) {
if(ta->a[k].occ == ((uint32_t)-1)) continue;
@@ -16771,7 +16861,7 @@ long long gap_fuzz, bub_label_t* b_mask_t)
asg_cleanup(sg);
// reduce_hamming_error(sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz);
reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, opt.ruIndex, NULL);
reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, opt.ruIndex, NULL, tipsLen);
output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex,
0.05, 0.9, max_hang, min_ovlp, gap_fuzz, 0, b_mask_t, NULL, NULL, NULL);
@@ -21282,6 +21372,45 @@ void output_hap_sc_graph(kvect_sec_t *ug, asg_t *sg, /**kvec_asg_arc_t_warp *arc
free(gfa_name);
}
scaf_res_t *gen_scaf_res_t(ug_opt_t *opt, ma_ug_t *sug, asg_t *rg, ma_ug_t *iug, ma_ug_t **res_ug)
{
ma_ug_t *ug = (iug?(iug):(ma_ug_gen(rg)));
scaf_res_t *sc = gen_contig_path(opt, rg, sug, ug);
if(res_ug) (*res_ug) = ug;
return sc;
}
void output_hap_sc_gpath(kvect_sec_t *ug, char* output_file_name, uint8_t flag, scaf_res_t *gpt, ma_ug_t *gpt_u, char *f_prefix)
{
char* gfa_name = (char*)malloc(strlen(output_file_name)+100);
sprintf(gfa_name, "%s.%s.p_ctg.path.paf", output_file_name, f_prefix?f_prefix:(flag==FATHER?"hap1":"hap2"));
FILE* output_file = fopen(gfa_name, "w");
fprintf(stderr, "Writing %s to disk... \n", gfa_name);
// ma_ug_seq(ug->ctg, sg, coverage_cut, sources, arcs, max_hang, min_ovlp, 0, 1);
// ma_ug_print(ug, sg, coverage_cut, sources, ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file);
ma_scg_gpath_print(ug, (flag==FATHER?"h1tg":"h2tg"), gpt, gpt_u, output_file);
fclose(output_file);
free(gfa_name);
}
void output_hap_bp_gpath(ug_opt_t *opt, ma_ug_t *ug, asg_t *sg, char* output_file_name, uint8_t flag, ma_ug_t *gpt_u, char *f_prefix)
{
scaf_res_t *gpt = gen_scaf_res_t(opt, ug, sg, gpt_u, NULL);
char* gfa_name = (char*)malloc(strlen(output_file_name)+100);
sprintf(gfa_name, "%s.%s.p_ctg.path.paf", output_file_name, f_prefix?f_prefix:(flag==FATHER?"hap1":"hap2"));
FILE* output_file = fopen(gfa_name, "w");
fprintf(stderr, "Writing %s to disk... \n", gfa_name);
ma_ug_gpath_print((flag==FATHER?"h1tg":"h2tg"), gpt, ug, gpt_u, output_file);
destroy_scaf_res_t(gpt);
fclose(output_file);
free(gfa_name);
}
void filter_set_kug(uint8_t* trio_flag, asg_t *rg, uint8_t *rf, kvec_asg_arc_t_warp *r_edges, float f_rate, ma_ug_t **ug)
{
@@ -21327,7 +21456,7 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long tipsLen, flo
long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang,
int min_ovlp, int is_bench, long long gap_fuzz, ug_opt_t *opt, bub_label_t* b_mask_t)
{
reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, opt->ruIndex, NULL);
reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, opt->ruIndex, NULL, tipsLen);
uint8_t *rf = NULL;
if(asm_opt.kpt_rate > 0) CALLOC(rf, sg->n_seq);
@@ -21421,7 +21550,7 @@ int64_t cal_exact_ug_o(dedup_idx_t *idx, ma_utg_t *u, uint64_t f)
{
if(u->n <= 0) return INT32_MIN;
uint64_t sv, *sa, sn, ev, *ea, en, si, ei, zn, *za, rev, k, nf = (uint64_t)-1, m, fn, nfn; ma_utg_t *z;
if(f == FATHER) nf = MOTHER; if(f == MOTHER) nf = FATHER;
if(f == FATHER) {nf = MOTHER;} if(f == MOTHER) {nf = FATHER;}
sv = u->a[0]>>32; sa = idx->ra + idx->ridx[sv>>1]; sn = idx->ridx[(sv>>1)+1]-idx->ridx[sv>>1];
ev = u->a[u->n-1]>>32; ea = idx->ra + idx->ridx[ev>>1]; en = idx->ridx[(ev>>1)+1] - idx->ridx[ev>>1];
for (si = 0; si < sn; si++) {
@@ -22705,8 +22834,8 @@ void gen_double_scaffold_gfa(ma_ug_t *ref, ma_ug_t *hu1, ma_ug_t *hu2, scaf_res_
r = u_trans_a(*os, k); rn = u_trans_n(*os, k); p = d = NULL; pn = dn = 0;
// fprintf(stderr, "[M::%s::] # k:%lu, rn::%lu\n", __func__, k, rn);
if(rn <= 1) continue;
for (z = 0, p = r; (z < rn) && (!(r[z].del)); z++); pn = z;
for (d = r + z; (z < rn) && (r[z].qs != ((uint32_t)-1)); z++); dn = z - pn;
for (z = 0, p = r; (z < rn) && (!(r[z].del)); z++){;} pn = z;
for (d = r + z; (z < rn) && (r[z].qs != ((uint32_t)-1)); z++){;} dn = z - pn;
// fprintf(stderr, "[M::%s::] # k:%lu, pn::%lu, dn::%lu\n", __func__, k, pn, dn);
if(pn <= 1) continue;///no scaf
s = e = (uint64_t)-1;
@@ -23241,6 +23370,17 @@ void merge_kvec_asg_arc_t_warp(kvec_asg_arc_t_warp *a, kvec_asg_arc_t_warp *b, k
memcpy(m->a.a + a->a.n, b->a.a, sizeof((*(b->a.a)))*b->a.n);
}
void output_noseq_gfa(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, ma_hit_t_alloc* sources, R_to_U* ruIndex)
{
char* gfa_name = (char*)malloc(strlen(output_file_name)+100);
sprintf(gfa_name, "%s.d_utg.noseq.gfa", output_file_name);
FILE* output_file = fopen(gfa_name, "w");
fprintf(stderr, "Writing %s to disk... \n", gfa_name);
ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file);
fclose(output_file);
free(gfa_name);
}
void output_trio_graph_joint(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long tipsLen, float tip_drop_ratio,
long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang,
@@ -23249,7 +23389,10 @@ int min_ovlp, long long gap_fuzz, bub_label_t* b_mask_t, ma_ug_t **rhu0, ma_ug_t
ma_ug_t *hu0 = NULL, *hu1 = NULL; kvec_asg_arc_t_warp arcs0, arcs1;
memset(&arcs0, 0, sizeof(arcs0)); memset(&arcs1, 0, sizeof(arcs1));
reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, ruIndex, NULL);
reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, ruIndex, NULL, tipsLen);
//debug
// print_debug_gfa(sg, NULL, opt->coverage_cut, "rd.hamm.debug", opt->sources, opt->ruIndex, opt->max_hang, opt->min_ovlp, 1, 0, 0);
// exit(1);
hu0 = output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources,
reverse_sources, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate,
@@ -23277,25 +23420,39 @@ int min_ovlp, long long gap_fuzz, bub_label_t* b_mask_t, ma_ug_t **rhu0, ma_ug_t
if((asm_opt.self_scaf) && (!rhu0) && (!rhu1)) {
kvect_sec_t *sc0 = NULL, *sc1 = NULL; ma_ug_t *mug = NULL; /**kv_u_trans_t *rs_trans = NULL;**/ kvec_asg_arc_t_warp ma; memset(&ma, 0, sizeof(ma));
scaf_res_t *gpt = NULL; ma_ug_t *gpt_u = NULL;
merge_kvec_asg_arc_t_warp(&arcs0, &arcs1, &ma); kv_destroy(arcs0.a); kv_destroy(arcs1.a);
gen_self_scaf(opt, hu0, hu1, sg, coverage_cut, sources, reverse_sources, ruIndex, tipsLen, tip_drop_ratio, stops_threshold, chimeric_rate, drop_ratio, max_hang, min_ovlp, b_mask_t, &sc0, &sc1, &mug, &ma, /**&rs_trans**/NULL);
ma_ug_destroy(hu0); ma_ug_destroy(hu1); kv_destroy(ma.a);
///output scaf trans
/**kv_destroy((*rs_trans)); kv_destroy(rs_trans->idx); free(rs_trans);**/
if(asm_opt.gpath) gpt = gen_scaf_res_t(opt, mug, sg, NULL, &gpt_u);
sc0->ctg = mug;
output_hap_sc_graph(sc0, sg, /**&arcs0,**/ coverage_cut, output_file_name, FATHER, sources, ruIndex, /**max_hang, min_ovlp,**/ NULL);
if(asm_opt.gpath) output_hap_sc_gpath(sc0, output_file_name, FATHER, gpt, gpt_u, NULL);
destory_kvect_sec_t(sc0); free(sc0);
sc1->ctg = mug;
output_hap_sc_graph(sc1, sg, /**&arcs1,**/ coverage_cut, output_file_name, MOTHER, sources, ruIndex, /**max_hang, min_ovlp,**/ NULL);
if(asm_opt.gpath) output_hap_sc_gpath(sc1, output_file_name, MOTHER, gpt, gpt_u, NULL);
destory_kvect_sec_t(sc1); free(sc1);
ma_ug_destroy(mug);
if(asm_opt.gpath) {
output_noseq_gfa(gpt_u, sg, coverage_cut, output_file_name, sources, ruIndex);
ma_ug_destroy(gpt_u);
}
ma_ug_destroy(mug); if(gpt) destroy_scaf_res_t(gpt);
} else {
ma_ug_t *gpt_u = NULL;
if((asm_opt.gpath) && ((!rhu0) || (!rhu1))) gpt_u = ma_ug_gen(sg);
if(!rhu0) {
output_hap_graph(hu0, sg, &arcs0, coverage_cut, output_file_name, FATHER, sources, ruIndex, max_hang, min_ovlp, NULL);
if(asm_opt.gpath) output_hap_bp_gpath(opt, hu0, sg, output_file_name, FATHER, gpt_u, NULL);
ma_ug_destroy(hu0);
} else {
(*rhu0) = hu0;
@@ -23305,11 +23462,17 @@ int min_ovlp, long long gap_fuzz, bub_label_t* b_mask_t, ma_ug_t **rhu0, ma_ug_t
if(!rhu1) {
output_hap_graph(hu1, sg, &arcs1, coverage_cut, output_file_name, MOTHER, sources, ruIndex, max_hang, min_ovlp, NULL);
if(asm_opt.gpath) output_hap_bp_gpath(opt, hu1, sg, output_file_name, MOTHER, gpt_u, NULL);
ma_ug_destroy(hu1);
} else {
(*rhu1) = hu1; // ma_ug_destroy(hu1);
}
kv_destroy(arcs1.a);
if(gpt_u) {
output_noseq_gfa(gpt_u, sg, coverage_cut, output_file_name, sources, ruIndex);
ma_ug_destroy(gpt_u);
}
}
}
@@ -23623,6 +23786,18 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s
return 0;
}
free(gfa_name);
if(asm_opt.write_pos_idx) {
int32_t m;
asm_opt.hom_cov_0 = asm_opt.hom_cov;
asm_opt.het_cov_0 = asm_opt.het_cov;
asm_opt.max_n_chain_0 = asm_opt.max_n_chain;
refresh_pt_idx(&ha_flt_tab, &ha_idx, NULL, &asm_opt, output_file_name, 0);
// load_pt_index(&ha_flt_tab, &ha_idx, NULL, &asm_opt, output_file_name);
m = asm_opt.hom_cov_0; asm_opt.hom_cov_0 = asm_opt.hom_cov; asm_opt.hom_cov = m;
m = asm_opt.het_cov_0; asm_opt.het_cov_0 = asm_opt.het_cov; asm_opt.het_cov = m;
m = asm_opt.max_n_chain_0; asm_opt.max_n_chain_0 = asm_opt.max_n_chain; asm_opt.max_n_chain = m;
}
return 1;
}
@@ -24693,6 +24868,17 @@ int asg_arc_del_trans_aux(asg_t *g, asg_t *aux, uint8_t *mark, int fuzz)
for (j = 0; j < nw1 && asg_arc_len(aw1[j]) + asg_arc_len(av0[i]) <= L; ++j)
if (mark[aw1[j].v]) mark[aw1[j].v] = 2;
}
for (i = 0; i < nv1; ++i) {
//w is an out-node of v
w = av1[i].v; if (mark[w] != 1) continue; ///w has already been reduced
nw0 = asg_arc_n(g, w); aw0 = asg_arc_a(g, w);
nw1 = asg_arc_n(aux, w); aw1 = asg_arc_a(aux, w);
for (j = 0; j < nw0 && asg_arc_len(aw0[j]) + asg_arc_len(av1[i]) <= L; ++j)
if (mark[aw0[j].v]) mark[aw0[j].v] = 2;
for (j = 0; j < nw1 && asg_arc_len(aw1[j]) + asg_arc_len(av1[i]) <= L; ++j)
if (mark[aw1[j].v]) mark[aw1[j].v] = 2;
}
//remove edges
for (i = 0; i < nv0; ++i) {
if (mark[av0[i].v] == 2) {
@@ -24700,6 +24886,7 @@ int asg_arc_del_trans_aux(asg_t *g, asg_t *aux, uint8_t *mark, int fuzz)
}
mark[av0[i].v] = 0;
}
for (i = 0; i < nv1; ++i) mark[av1[i].v] = 0;
}
if (n_reduced) {
@@ -25411,15 +25598,272 @@ void rd_hamming_symm_simple(rd_hamming_t *s, uint32_t st, uint32_t ed) // callba
}
}
uint32_t rd_hamming_symm_simple0(buf_t *b, asg_t *ref, asg_t *g, uint32_t st, uint32_t ed, uint64_t max_dist, uint64_t *r_max_dist) // callback for kt_for()
uint32_t pp_cut(uint32_t vi0, uint32_t v0, asg_t *g, asg_t *ref, asg32_v *sa, asg32_v *sb, uint32_t max_ext)
{
//e078bbb4-47b8-49af-9221-889e67f7097a 0 47735 id:i:1477089 (w)
//f552c01e-a9ea-47f2-b1ab-a56f6afa3121 0 42133 id:i:1477068 (v)
// if((vi0>>1) == 1476811 || (vi0>>1) == 1479120) {
// fprintf(stderr, "\n\ns1s[M::%s] S\t%.*s(%c)\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, v0>>1), Get_NAME(R_INF, v0>>1), "+-"[v0&1]);
// }
uint32_t k, i, nv, nw, kv, kw, w, v = v0, ntip = 0, f; asg_arc_t *av, *aw;
nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
for (k = 0; k < nv && av[k].del; ++k);
if(k >= nv) return 1;///no new arcs
// if((vi0>>1) == 1476811 || (vi0>>1) == 1479120) {
// fprintf(stderr, "s2s[M::%s] S\t%.*s(%c)\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, v0>>1), Get_NAME(R_INF, v0>>1), "+-"[v0&1]);
// }
sa->n = sb->n = 0;
nv = asg_arc_n(ref, v); av = asg_arc_a(ref, v);
for (k = kv = 0; k < nv; ++k) {
if(av[k].del) continue;
kv++; w = av[k].v^1;
kv_push(u_int32_t, *sb, av-ref->arc + k + g->n_arc);
av[k].del = 1;
nw = asg_arc_n(ref, w); aw = asg_arc_a(ref, w);
for (i = kw = f = 0; i < nw; i++) {
if(aw[i].del) continue;
if(aw[i].v == (v^1)) {
kv_push(u_int32_t, *sb, aw-ref->arc + i + g->n_arc);
aw[i].del = 1; f = 1;
} else {
// if(((vi0>>1) == 1476811 || (vi0>>1) == 1479120) && ((v>>1) == 1477068) && ((w>>1) == 1477089)) {
// fprintf(stderr, "arc_0[M::%s] L\t%.*s(%c)\t%.*s(%c)\tkw::%u\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, aw[i].ul>>33), Get_NAME(R_INF, aw[i].ul>>33), "+-"[(aw[i].ul>>32)&1],
// (int)Get_NAME_LENGTH(R_INF, aw[i].v>>1), Get_NAME(R_INF, aw[i].v>>1), "+-"[aw[i].v&1], kw);
// }
kw++;
}
}
assert(f);
nw = asg_arc_n(g, w); aw = asg_arc_a(g, w);
for (i = 0; i < nw; i++) {
if(aw[i].del) continue;
// if(((vi0>>1) == 1476811 || (vi0>>1) == 1479120) && ((v>>1) == 1477068) && ((w>>1) == 1477089)) {
// fprintf(stderr, "arc_1[M::%s] L\t%.*s(%c)\t%.*s(%c)\tkw::%u\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, aw[i].ul>>33), Get_NAME(R_INF, aw[i].ul>>33), "+-"[(aw[i].ul>>32)&1],
// (int)Get_NAME_LENGTH(R_INF, aw[i].v>>1), Get_NAME(R_INF, aw[i].v>>1), "+-"[aw[i].v&1], kw);
// }
kw++;
}
// if((vi0>>1) == 1476811 || (vi0>>1) == 1479120) {
// fprintf(stderr, "init_a[M::%s] L\t%.*s(%c)\t%.*s(%c)\tkw::%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], kw);
// }
if(!kw) {///a new tip
// if((vi0>>1) == 1476811 || (vi0>>1) == 1479120) {
// fprintf(stderr, "[M::%s] L\t%.*s(%c)\t%.*s(%c)\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, v0>>1), Get_NAME(R_INF, v0>>1), "+-"[v0&1],
// (int)Get_NAME_LENGTH(R_INF, w>>1), Get_NAME(R_INF, w>>1), "+-"[w&1]);
// }
kv_push(uint32_t, *sa, (w^1)); ntip++;
}
}
while (sa->n) {
v = kv_pop(*sa);
nv = asg_arc_n(ref, v); av = asg_arc_a(ref, v);
for (k = kv = 0; k < nv; ++k) {
if(av[k].del) continue;
kv_push(u_int32_t, *sb, av-ref->arc + k + g->n_arc);
av[k].del = 1; kv++; w = av[k].v^1; f = 0;
nw = asg_arc_n(ref, w); aw = asg_arc_a(ref, w);
for (i = kw = 0; i < nw; i++) {
if(aw[i].del) continue;
if(aw[i].v == (v^1)) {
kv_push(u_int32_t, *sb, aw-ref->arc + i + g->n_arc);
aw[i].del = 1; f = 1;
} else {
// if(((vi0>>1) == 1476811 || (vi0>>1) == 1479120) && ((v>>1) == 1477068) && ((w>>1) == 1477089)) {
// fprintf(stderr, "arc_0_a[M::%s] L\t%.*s(%c)\t%.*s(%c)\tkw::%u\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, aw[i].ul>>33), Get_NAME(R_INF, aw[i].ul>>33), "+-"[(aw[i].ul>>32)&1],
// (int)Get_NAME_LENGTH(R_INF, aw[i].v>>1), Get_NAME(R_INF, aw[i].v>>1), "+-"[aw[i].v&1], kw);
// }
kw++;
}
}
assert(f);
nw = asg_arc_n(g, w); aw = asg_arc_a(g, w);
for (i = 0; i < nw; i++) {
if(aw[i].del) continue;
// if(((vi0>>1) == 1476811 || (vi0>>1) == 1479120) && ((v>>1) == 1477068) && ((w>>1) == 1477089)) {
// fprintf(stderr, "arc_1_a[M::%s] L\t%.*s(%c)\t%.*s(%c)\tkw::%u\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, aw[i].ul>>33), Get_NAME(R_INF, aw[i].ul>>33), "+-"[(aw[i].ul>>32)&1],
// (int)Get_NAME_LENGTH(R_INF, aw[i].v>>1), Get_NAME(R_INF, aw[i].v>>1), "+-"[aw[i].v&1], kw);
// }
kw++;
}
// if((vi0>>1) == 1476811 || (vi0>>1) == 1479120) {
// fprintf(stderr, "init_b[M::%s] L\t%.*s(%c)\t%.*s(%c)\tkw::%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], kw);
// }
if(!kw) {
// if((vi0>>1) == 1476811 || (vi0>>1) == 1479120) {
// fprintf(stderr, "[M::%s] L\t%.*s(%c)\t%.*s(%c)\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, v0>>1), Get_NAME(R_INF, v0>>1), "+-"[v0&1],
// (int)Get_NAME_LENGTH(R_INF, w>>1), Get_NAME(R_INF, w>>1), "+-"[w&1]);
// }
kv_push(uint32_t, *sa, (w^1)); ntip++;
}
}
nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
for (k = 0; k < nv; ++k) {
if(av[k].del) continue;
kv_push(u_int32_t, *sb, av-g->arc + k);
av[k].del = 1; kv++; w = av[k].v^1; f = 0;
nw = asg_arc_n(ref, w); aw = asg_arc_a(ref, w);
for (i = kw = 0; i < nw; i++) {
if(aw[i].del) continue;
// if(((vi0>>1) == 1476811 || (vi0>>1) == 1479120) && ((v>>1) == 1477068) && ((w>>1) == 1477089)) {
// fprintf(stderr, "arc_2_a[M::%s] L\t%.*s(%c)\t%.*s(%c)\tkw::%u\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, aw[i].ul>>33), Get_NAME(R_INF, aw[i].ul>>33), "+-"[(aw[i].ul>>32)&1],
// (int)Get_NAME_LENGTH(R_INF, aw[i].v>>1), Get_NAME(R_INF, aw[i].v>>1), "+-"[aw[i].v&1], kw);
// }
kw++;
}
nw = asg_arc_n(g, w); aw = asg_arc_a(g, w);
for (i = 0; i < nw; i++) {
if(aw[i].del) continue;
if(aw[i].v == (v^1)) {
kv_push(u_int32_t, *sb, aw-g->arc + i);
aw[i].del = 1; f = 1;
} else {
// if(((vi0>>1) == 1476811 || (vi0>>1) == 1479120) && ((v>>1) == 1477068) && ((w>>1) == 1477089)) {
// fprintf(stderr, "arc_2_b[M::%s] L\t%.*s(%c)\t%.*s(%c)\tkw::%u\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, aw[i].ul>>33), Get_NAME(R_INF, aw[i].ul>>33), "+-"[(aw[i].ul>>32)&1],
// (int)Get_NAME_LENGTH(R_INF, aw[i].v>>1), Get_NAME(R_INF, aw[i].v>>1), "+-"[aw[i].v&1], kw);
// }
kw++;
}
}
assert(f);
// if((vi0>>1) == 1476811 || (vi0>>1) == 1479120) {
// fprintf(stderr, "init_c[M::%s] L\t%.*s(%c)\t%.*s(%c)\tkw::%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], kw);
// }
if(!kw) {
// if((vi0>>1) == 1476811 || (vi0>>1) == 1479120) {
// fprintf(stderr, "[M::%s] L\t%.*s(%c)\t%.*s(%c)\n", __func__,
// (int)Get_NAME_LENGTH(R_INF, v0>>1), Get_NAME(R_INF, v0>>1), "+-"[v0&1],
// (int)Get_NAME_LENGTH(R_INF, w>>1), Get_NAME(R_INF, w>>1), "+-"[w&1]);
// }
kv_push(uint32_t, *sa, (w^1)); ntip++;
}
}
}
for (i = 0; i < sb->n; i++) {
if(sb->a[i] < g->n_arc) {
assert(g->arc[sb->a[i]].del);
g->arc[sb->a[i]].del = 0;
} else {
assert(ref->arc[sb->a[i] - g->n_arc].del);
ref->arc[sb->a[i] - g->n_arc].del = 0;
}
}
if(ntip > 0 && ntip >= max_ext) return 0;
return 1;
}
uint32_t refine_rd_hamming_symm(asg_t *g, asg_t *ref, uint32_t v0, uint32_t v1, buf_t *z, uint32_t max_ext)
{
// if((v0>>1) == 1476811 || (v0>>1) == 1479120) {
// fprintf(stderr, "[M::%s]\n", __func__);
// }
uint32_t i, ncut = 0;
uint32_t v, w, nv, k; asg_arc_t *av;
if (g->seq[v0>>1].del) return 0; // already deleted
asg32_v ai, bi, ci;
copy_asg_arr(ai, z->S); copy_asg_arr(bi, z->T); copy_asg_arr(ci, z->b);
ai.n = bi.n = ci.n = 0;
kv_push(uint32_t, ai, v0);
while(ai.n) {
v = kv_pop(ai);
if(z->a[v].s) continue;
z->a[v].s = 1;
kv_push(uint32_t, bi, v); // save it for revert
nv = asg_arc_n(ref, v); av = asg_arc_a(ref, v);
for (k = 0; k < nv; ++k) { // loop through v's neighbors
if (av[k].del) continue;
w = av[k].v;
if(z->a[w].s || w == v1) continue;
kv_push(uint32_t, ai, w);
}
}
ncut = 1;
while (ncut) {
for (i = ncut = 0; i < bi.n; ++i) { // clear the states of visited vertices
v = bi.a[i]; z->a[v].s = 0;
if(v == v0 || v == v1) continue;
if(!pp_cut(v0, v, g, ref, &ai, &ci, max_ext)) {
nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
for (k = 0; k < nv; ++k) { // loop through v's neighbors
if (av[k].del) continue;
w = av[k].v; av[k].del = 1;
asg_arc_del(g, av[k].v^1, (av[k].ul>>32)^1, 1); ncut++;
// fprintf(stderr, "L\t%.*s(%c)\t%.*s(%c)\n",
// (int)Get_NAME_LENGTH(R_INF, (av[k].ul>>33)), Get_NAME(R_INF, (av[k].ul>>33)), "+-"[(av[k].ul>>32)&1],
// (int)Get_NAME_LENGTH(R_INF, (av[k].v>>1)), Get_NAME(R_INF, (av[k].v>>1)), "+-"[av[k].v&1]);
}
}
if(!pp_cut(v0, v^1, g, ref, &ai, &ci, max_ext)) {
nv = asg_arc_n(g, v^1); av = asg_arc_a(g, v^1);
for (k = 0; k < nv; ++k) { // loop through v's neighbors
if (av[k].del) continue;
w = av[k].v; av[k].del = 1;
asg_arc_del(g, av[k].v^1, (av[k].ul>>32)^1, 1); ncut++;
// fprintf(stderr, "L\t%.*s(%c)\t%.*s(%c)\n",
// (int)Get_NAME_LENGTH(R_INF, (av[k].ul>>33)), Get_NAME(R_INF, (av[k].ul>>33)), "+-"[(av[k].ul>>32)&1],
// (int)Get_NAME_LENGTH(R_INF, (av[k].v>>1)), Get_NAME(R_INF, (av[k].v>>1)), "+-"[av[k].v&1]);
}
}
}
}
ai.n = bi.n = ci.n = 0;
copy_asg_arr(z->S, ai); copy_asg_arr(z->T, bi); copy_asg_arr(z->b, ci);
return ncut;
}
uint32_t rd_hamming_symm_simple0(buf_t *b, asg_t *ref, asg_t *g, uint32_t st, uint32_t ed, uint64_t max_dist, uint64_t *r_max_dist, uint32_t max_ext) // callback for kt_for()
{
// if((st>>1) == 1476811 || (st>>1) == 1479120) {
// fprintf(stderr, "[M::%s]\tst>>1::%u(%c)\ted>>1::%u(%c)\n", __func__, st>>1, "+-"[st&1], ed>>1, "+-"[ed&1]);
// }
double step = 0.2, cuttoff;
uint32_t p, k, ncut;
p = rd_hm_bub(g, ref, st, max_dist, b);
if(p) {
// if((st>>1) == 1476811 || (st>>1) == 1479120) {
// fprintf(stderr, "-a-[M::%s]\tst>>1::%u(%c)\ted>>1::%u(%c)\n", __func__, st>>1, "+-"[st&1], ed>>1, "+-"[ed&1]);
// }
assert(b->S.a[0] == (ed^1));
if(r_max_dist) (*r_max_dist) = max_dist;
refine_rd_hamming_symm(g, ref, st, ed^1, b, max_ext);
return 1;
}
///recalculate max_dist
@@ -25437,8 +25881,12 @@ uint32_t rd_hamming_symm_simple0(buf_t *b, asg_t *ref, asg_t *g, uint32_t st, ui
max_dist += ref->seq[b->S.a[0]>>1].len;
p = rd_hm_bub(g, ref, st, max_dist, b);
if(p) {
// if((st>>1) == 1476811 || (st>>1) == 1479120) {
// fprintf(stderr, "-b-[M::%s]\tst>>1::%u(%c)\ted>>1::%u(%c)\n", __func__, st>>1, "+-"[st&1], ed>>1, "+-"[ed&1]);
// }
assert(b->S.a[0] == (ed^1));
if(r_max_dist) (*r_max_dist) = max_dist;
refine_rd_hamming_symm(g, ref, st, ed^1, b, max_ext);
return 1;
}
@@ -25447,8 +25895,12 @@ uint32_t rd_hamming_symm_simple0(buf_t *b, asg_t *ref, asg_t *g, uint32_t st, ui
ncut = rd_hm_drop(g, ref, st, ed^1, cuttoff, 1, b);
p = rd_hm_bub(g, ref, st, max_dist, b);
if(p) {
// if((st>>1) == 1476811 || (st>>1) == 1479120) {
// fprintf(stderr, "-c-[M::%s]\tst>>1::%u(%c)\ted>>1::%u(%c)\n", __func__, st>>1, "+-"[st&1], ed>>1, "+-"[ed&1]);
// }
assert(b->S.a[0] == (ed^1));
if(r_max_dist) (*r_max_dist) = max_dist;
refine_rd_hamming_symm(g, ref, st, ed^1, b, max_ext);
return 1;
}
@@ -25456,8 +25908,12 @@ uint32_t rd_hamming_symm_simple0(buf_t *b, asg_t *ref, asg_t *g, uint32_t st, ui
ncut = rd_hm_drop(g, ref, st, ed^1, cuttoff, 0, b);
p = rd_hm_bub(g, ref, st, max_dist, b);
if(p) {
// if((st>>1) == 1476811 || (st>>1) == 1479120) {
// fprintf(stderr, "-d-[M::%s]\tst>>1::%u(%c)\ted>>1::%u(%c)\n", __func__, st>>1, "+-"[st&1], ed>>1, "+-"[ed&1]);
// }
assert(b->S.a[0] == (ed^1));
if(r_max_dist) (*r_max_dist) = max_dist;
refine_rd_hamming_symm(g, ref, st, ed^1, b, max_ext);
return 1;
}
if(!ncut) break;
@@ -25470,8 +25926,55 @@ uint32_t rd_hamming_symm_simple0(buf_t *b, asg_t *ref, asg_t *g, uint32_t st, ui
return 0;
}
uint32_t prt_arc_status(asg_t *g, uint32_t v0, uint32_t w0, const char *cmd)
{
uint32_t v, w, i, nv; asg_arc_t *av;
v = v0<<1; w = w0;
nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
for (i = 0; i < nv; ++i) {
if (((av[i].v>>1) == w) && (!av[i].del)) {
fprintf(stderr, "%s\tL\t%.*s(%c)\t%.*s(%c)\n", cmd,
(int)Get_NAME_LENGTH(R_INF, (av[i].ul>>33)), Get_NAME(R_INF, (av[i].ul>>33)), "+-"[(av[i].ul>>32)&1],
(int)Get_NAME_LENGTH(R_INF, (av[i].v>>1)), Get_NAME(R_INF, (av[i].v>>1)), "+-"[av[i].v&1]);
}
}
v = (v0<<1) + 1; w = w0;
nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
for (i = 0; i < nv; ++i) {
if (((av[i].v>>1) == w) && (!av[i].del)) {
fprintf(stderr, "%s\tL\t%.*s(%c)\t%.*s(%c)\n", cmd,
(int)Get_NAME_LENGTH(R_INF, (av[i].ul>>33)), Get_NAME(R_INF, (av[i].ul>>33)), "+-"[(av[i].ul>>32)&1],
(int)Get_NAME_LENGTH(R_INF, (av[i].v>>1)), Get_NAME(R_INF, (av[i].v>>1)), "+-"[av[i].v&1]);
}
}
v = w0<<1; w = v0;
nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
for (i = 0; i < nv; ++i) {
if (((av[i].v>>1) == w) && (!av[i].del)) {
fprintf(stderr, "%s\tL\t%.*s(%c)\t%.*s(%c)\n", cmd,
(int)Get_NAME_LENGTH(R_INF, (av[i].ul>>33)), Get_NAME(R_INF, (av[i].ul>>33)), "+-"[(av[i].ul>>32)&1],
(int)Get_NAME_LENGTH(R_INF, (av[i].v>>1)), Get_NAME(R_INF, (av[i].v>>1)), "+-"[av[i].v&1]);
}
}
v = (w0<<1) + 1; w = v0;
nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
for (i = 0; i < nv; ++i) {
if (((av[i].v>>1) == w) && (!av[i].del)) {
fprintf(stderr, "%s\tL\t%.*s(%c)\t%.*s(%c)\n", cmd,
(int)Get_NAME_LENGTH(R_INF, (av[i].ul>>33)), Get_NAME(R_INF, (av[i].ul>>33)), "+-"[(av[i].ul>>32)&1],
(int)Get_NAME_LENGTH(R_INF, (av[i].v>>1)), Get_NAME(R_INF, (av[i].v>>1)), "+-"[av[i].v&1]);
}
}
return 1;
}
void reduce_hamming_error_adv(ma_ug_t *iug, asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub)
int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub, uint32_t max_ext)
{
double index_time = yak_realtime();
ma_ug_t *ug = NULL; ug = (iug)?(iug):(ma_ug_gen_primary(sg, PRIMARY_LABLE));
@@ -25534,13 +26037,23 @@ int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub)
free(bs_flag); if(!iug) ma_ug_destroy(ug);
if(sv.n > 0) {
ig->n_seq = ig->m_seq = sg->n_seq;
MALLOC(ig->seq, ig->n_seq);
memcpy(ig->seq, sg->seq, (sizeof((*(ig->seq)))*ig->n_seq));
asg_cleanup(ig); asg_arc_del_trans_aux(ig, sg, vis_flag, gap_fuzz);
REALLOC(b.a, (ig->n_seq<<1)); memset(b.a, 0, sizeof((*(b.a)))*(ig->n_seq<<1));
for (i = 0; i < sv.n; i++) rd_hamming_symm_simple0(&b, sg, ig, sv.a[i]>>32, (uint32_t)(sv.a[i]), max_dist, NULL);
// prt_arc_status(ig, 1477050, 1477089, "init");
for (i = 0; i < sv.n; i++) {
rd_hamming_symm_simple0(&b, sg, ig, sv.a[i]>>32, (uint32_t)(sv.a[i]), max_dist, NULL, max_ext);
// fprintf(stderr, "[M::%s]\tst>>1::%lu(%c)\ted>>1::%u(%c)\n", __func__, (sv.a[i]>>32)>>1, "+-"[(sv.a[i]>>32)&1], ((uint32_t)(sv.a[i]))>>1, "+-"[(uint32_t)(sv.a[i])&1]);
// prt_arc_status(ig, 1477050, 1477089, "mm");
}
// prt_arc_status(ig, 1477050, 1477089, "end");
/**
rd_hamming_t aux_t; memset((&aux_t), 0, sizeof(aux_t));
aux_t.n_thread = 1; // aux_t.n_thread = asm_opt.thread_num;
@@ -25561,6 +26074,8 @@ int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub)
}
free(sv.a); free(vis_flag);
// prt_arc_status(ig, 1477050, 1477089, "fin");
for (i = n_pop = 0; i < ig->n_arc; i++) {
if(ig->arc[i].del) continue;
p = asg_arc_pushp(sg); *p = (ig->arc[i]); n_pop++;
@@ -25570,9 +26085,15 @@ int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub)
sg->idx = 0;
sg->is_srt = 0;
asg_cleanup(sg);
// prt_arc_status(sg, 1477050, 1477089, "ff0");
asg_symm(sg);
// prt_arc_status(sg, 1477050, 1477089, "ff1");
asg_arc_del_trans(sg, gap_fuzz);
// prt_arc_status(sg, 1477050, 1477089, "ff2");
}
// prt_arc_status(sg, 1477050, 1477089, "ffe");
free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); asg_destroy(ig);
fprintf(stderr, "[M::%s::%.3f] # inserted edges: %u, # fixed bubbles: %u\n",
__func__, yak_realtime() - index_time, sg->n_arc - n_arc_0, fix_bub);
@@ -25819,8 +26340,7 @@ hap_cov_t *cov, uint32_t is_update_chain, uint32_t keep_d, utg_trans_t *o)
///v is a node that all incoming edges have been visited
///d is the distance from v0 to v
uint32_t v = kv_pop(b->S), d = b->a[v].d, c = b->a[v].c, m = b->a[v].m, nc = b->a[v].nc, np = b->a[v].np;
uint32_t nv = asg_arc_n(g, v);
asg_arc_t *av = asg_arc_a(g, v);
uint32_t nv = asg_arc_n(g, v); asg_arc_t *av = asg_arc_a(g, v);
///why we have this assert?
///assert(nv > 0);
///all out-edges of v
@@ -26530,7 +27050,7 @@ static inline void asg_arc_rest(asg_t *des, asg_t *src, uint32_t v0, uint32_t w0
assert(nv);
for (i = 0; i < nv; ++i) assert(av[i].del);
for (i = 0; i < nv && av[i].ul < arc->ul; ++i);
if(i >= nv) i = nv - 1; av[i] = *arc;
if(i >= nv) {i = nv - 1;} av[i] = *arc;
v = w0^1; w = v0^1;
av = asg_arc_a(src, v); nv = asg_arc_n(src, v);
@@ -26544,7 +27064,7 @@ static inline void asg_arc_rest(asg_t *des, asg_t *src, uint32_t v0, uint32_t w0
assert(nv);
for (i = 0; i < nv; ++i) assert(av[i].del);
for (i = 0; i < nv && av[i].ul < arc->ul; ++i);
if(i >= nv) i = nv - 1; av[i] = *arc;
if(i >= nv) {i = nv - 1;} av[i] = *arc;
rv = ((v0&1)?(Uc_beg(ug->u.a[v0>>1])^1):(Uc_end(ug->u.a[v0>>1])^1));
@@ -26711,7 +27231,7 @@ uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov, utg_trans_t *o,
}
is_update = rd_hamming_symm_simple0(b, ug->g, fg->g, v0, v1^1, max_dist, &max_dist);
is_update = rd_hamming_symm_simple0(b, ug->g, fg->g, v0, v1^1, max_dist, &max_dist, (asm_opt.max_short_tip*2));
// if(((v0>>1) == 7536) && ((v1>>1) == 99223)) {
// fprintf(stderr, "v0::utg%.6lul(%c)\tv1::utg%.6lul(%c)\n", (v0>>1) + 1, "+-"[(v0&1)], (v1>>1) + 1, "+-"[(v1&1)]);
@@ -39241,7 +39761,7 @@ ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t,
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);
// prt_dbg_rid_ovlp(src, -1, (char*)"c85c2e91-0490-438b-977b-b7d056973996", "7");
// 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);
@@ -39265,6 +39785,7 @@ ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t,
init_bub_label_t(b_mask_t, MIN(10, asm_opt.thread_num), sg->n_seq);
asm_opt.coverage = get_coverage(src, *cov, n_read);
// prt_specfic_sge(sg, 22708, 22646, "--#--");
return sg;
}
@@ -39611,8 +40132,12 @@ ma_sub_t **coverage_cut_ptr, uint8_t *cmk, int debug_g)
write_debug_graph(sg, sources, coverage_cut, output_file_name, n_read, reverse_sources, ruIndex);
debug_gfa:;
}**/
if(asm_opt.fn_bin_poy)
if(asm_opt.fn_chr_bin) {
if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION;
output_chr_bin_trio(sg, coverage_cut, o_file, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex,
0.05, 0.9, max_hang_length, mini_overlap_length, gap_fuzz, 0, &b_mask_t);
}
else if(asm_opt.fn_bin_poy)
{
if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION;
output_poly_trio(sg, coverage_cut, o_file, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex,
Process_Read.cpp
+40 -21
View File
@@ -46,14 +46,20 @@ void init_All_reads(All_reads* r)
void destory_All_reads(All_reads* r)
{
uint64_t i = 0;
for (i = 0; i < r->total_reads; i++) {
for (i = 0; i < r->tqn; i++) {
if (r->N_site[i]) free(r->N_site[i]);
if (r->read_sperate[i]) free(r->read_sperate[i]);
if (r->paf && r->paf[i].buffer) free(r->paf[i].buffer);
if (r->reverse_paf && r->reverse_paf[i].buffer) free(r->reverse_paf[i].buffer);
if(r->rsc && r->rsc[i]) free(r->rsc[i]);
///if (r->pb_regions) kv_destroy(r->pb_regions[i].a);
}
for (; i < r->total_reads; i++) {
if (r->N_site[i]) free(r->N_site[i]);
if (r->read_sperate[i]) free(r->read_sperate[i]);
if (r->paf && r->paf[i].buffer) free(r->paf[i].buffer);
if (r->reverse_paf && r->reverse_paf[i].buffer) free(r->reverse_paf[i].buffer);
}
free(r->paf);
free(r->reverse_paf);
free(r->N_site);
@@ -110,10 +116,11 @@ void write_All_reads(All_reads* r, char* read_file_name)
fwrite(&(asm_opt.hom_cov), sizeof(asm_opt.hom_cov), 1, fp);
fwrite(&(asm_opt.het_cov), sizeof(asm_opt.het_cov), 1, fp);
uint64_t mm = 1;
uint64_t mm = 2;///1;
if(asm_opt.is_sc) {
fwrite(&mm, sizeof(mm), 1, fp);
for (i = 0; i < r->total_reads; i++) {
fwrite(&(r->tqn), sizeof(r->tqn), 1, fp);
for (i = 0; i < r->tqn; i++) {
fwrite(r->rsc[i], sizeof(uint8_t), ((r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0)), fp);
}
}
@@ -133,6 +140,7 @@ int load_All_reads(All_reads* r, char* read_file_name)
free(index_name);
return 0;
}
// fprintf(stderr, "[M::%s]\tindex_name::%s\n", __func__, index_name);
int local_adapterLen;
int f_flag;
f_flag = fread(&local_adapterLen, sizeof(local_adapterLen), 1, fp);
@@ -213,9 +221,17 @@ int load_All_reads(All_reads* r, char* read_file_name)
uint64_t mm = 0;
if (!feof(fp)) {
if((fread(&mm, sizeof(mm), 1, fp)) && (mm == 1)) {
MALLOC(r->rsc, r->total_reads);
for (i = 0; i < r->total_reads; i++) {
if((fread(&mm, sizeof(mm), 1, fp)) && (mm == 1 || mm == 2)) {
if(mm == 1) {
mm = r->total_reads;
} else {
assert(mm == 2);
fread(&mm, sizeof(mm), 1, fp);
}
r->tqn = mm;
MALLOC(r->rsc, r->tqn);
for (i = 0; i < r->tqn; i++) {
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);
}
@@ -274,8 +290,8 @@ int append_All_reads(All_reads* r, char *idx, uint32_t id)
int local_adapterLen, f_flag;
f_flag = fread(&local_adapterLen, sizeof(local_adapterLen), 1, fp);
if(local_adapterLen != asm_opt.adapterLen) {
fprintf(stderr, "the adapterLen of index is: %d, but the adapterLen set by user is: %d\n",
local_adapterLen, asm_opt.adapterLen);
fprintf(stderr, "[M::%s] the adapterLen of index is: %d, but the adapterLen set by user is: %d\n",
__func__, local_adapterLen, asm_opt.adapterLen);
exit(1);
}
uint64_t index_size0, name_index_size0, total_reads0, total_reads_bases0, total_name_length0;
@@ -435,13 +451,18 @@ void malloc_All_reads(All_reads* r)
memcpy(r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads);
r->read_sperate = (uint8_t**)malloc(sizeof(uint8_t*)*r->total_reads);
if(asm_opt.is_sc) MALLOC(r->rsc, r->total_reads);
if(asm_opt.is_sc && r->tqn) MALLOC(r->rsc, r->tqn);
assert(r->tqn <= r->total_reads);
long long i = 0;
for (i = 0; i < (long long)r->total_reads; i++)
{
r->read_sperate[i] = (uint8_t*)malloc(sizeof(uint8_t)*(r->read_length[i]/4+1));
if(r->rsc) MALLOC(r->rsc[i], (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0));
uint64_t i = 0;
if(r->rsc) {
for (i = 0; i < r->tqn; i++) {
MALLOC(r->read_sperate[i], (r->read_length[i]/4+1));
MALLOC(r->rsc[i], (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0));
}
}
for (; i < r->total_reads; i++) {
MALLOC(r->read_sperate[i], (r->read_length[i]/4+1));
}
r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
@@ -450,8 +471,7 @@ void malloc_All_reads(All_reads* r)
r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
///r->pb_regions = (kvec_t_u64_warp*)malloc(r->total_reads*sizeof(kvec_t_u64_warp));
for (i = 0; i < (long long)r->total_reads; i++)
{
for (i = 0; i < r->total_reads; i++) {
r->second_round_cigar[i].size = r->cigars[i].size = 0;
r->second_round_cigar[i].length = r->cigars[i].length = 0;
r->second_round_cigar[i].record = r->cigars[i].record = NULL;
@@ -915,14 +935,14 @@ void ha_compress_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn, u
int64_t retrive_bqual(asg8_v *dv, uint8_t *ds, uint64_t id, int64_t s, int64_t e, uint8_t rev, int64_t bitn)
{
int64_t rl = Get_READ_LENGTH(R_INF, id), l;
if(s < 0) s = 0; if(e < 0) e = rl;
if(s < 0) {s = 0;} if(e < 0) {e = rl;}
if(s >= e || e > rl) return -1;
uint8_t *da = NULL, *src = Get_QUAL(R_INF, id), mm = (((uint8_t)1)<<bitn)-1, mlf = 8 - bitn, mrf;
int64_t bitr = 8/bitn, dk, sk, swk;
l = e - s;
if(dv) {
kv_resize(uint8_t, *dv, ((uint64_t)l)); da = dv->a;
kv_resize(uint8_t, *dv, ((uint64_t)l)); da = dv->a; dv->n = l;
} else {
da = ds;
}
@@ -969,7 +989,6 @@ int64_t retrive_bqual(asg8_v *dv, uint8_t *ds, uint64_t id, int64_t s, int64_t e
assert(dk == 0);
}
dv->n = l;
return l;
}
@@ -1499,7 +1518,7 @@ void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str,
p->rlen = str_l;
// fprintf(stderr, "str_l->%ld, str->%u\n", str_l, str?1:0);
if(o == NULL || on == 0) on = 0; en = 0;
if(o == NULL || on == 0) {on = 0;} en = 0;
for (i = on-1, st = et = str_l; i >= 0; i--) {
z = &(o[i]);
if(z->el) {
Process_Read.h
+2
View File
@@ -133,8 +133,10 @@ typedef struct
uint64_t* name_index;
uint64_t name_index_size;
uint64_t total_reads;
uint64_t tqn;
uint64_t total_reads_bases;
uint64_t total_name_length;
uint64_t tr[2];
Compressed_Cigar_record* cigars;
Compressed_Cigar_record* second_round_cigar;
Trio.cpp
+93
View File
@@ -388,6 +388,31 @@ inline void phrase_hstatus(char *s, char **rname, uint32_t *hid)
*hid = atoi(id);
}
inline void phrase_hchar(char *s, char **rname, uint32_t *hid)
{
*rname = NULL; *hid = (uint32_t)-1;
uint32_t tot = 0, k, l, z, sl = strlen(s);
for (k = 1, l = 0; k <= sl; k++) {
if((k == sl) || (s[k] == '\t') || (s[k] == ' ')) {
if(k > l) {
s[k] = 0;
if(tot == 0) {
*rname = s + l;
} else if(tot == 1) {
for (z = l; (z < k) && (s[z] >= '0') && (s[z] <= '9'); ++z);
if(z < k) {
fprintf(stderr, "ERROR: wrong hap id\n");
return;
}
*hid = atoi(s + l);
}
tot++;
}
l = k + 1;
}
}
}
uint32_t *ha_polybin_list(const hifiasm_opt_t *opt)
{
int64_t i;
@@ -447,6 +472,74 @@ uint32_t *ha_polybin_list(const hifiasm_opt_t *opt)
return ss;
}
uint32_t *ha_charbin_list(const hifiasm_opt_t *opt, uint8_t **idx, uint32_t *idx_n)
{
int64_t i;
khint_t k;
cstr_ht_t *h; (*idx) = NULL; *idx_n = 0;
assert(R_INF.total_reads < (uint32_t)-1);
h = cstr_ht_init();
for (i = 0; i < (int64_t)R_INF.total_reads; ++i) {
int absent;
char *str = (char*)calloc(Get_NAME_LENGTH(R_INF, i) + 1, 1);
strncpy(str, Get_NAME(R_INF, i), Get_NAME_LENGTH(R_INF, i));
k = cstr_ht_put(h, str, &absent);
if (absent) kh_val(h, k) = i;
}
fprintf(stderr, "[M::%s::%.3f*%.2f] created the hash table for read names\n", __func__, yak_realtime(), yak_cpu_usage());
gzFile fp;
kstream_t *ks;
kstring_t str = {0,0,0};
char *rname = NULL;
uint32_t hid, mhid = 0, *ss = NULL;
int dret;
int64_t n_tot = 0, n_bin = 0;
fp = gzopen(opt->fn_chr_bin, "r");
if (fp == 0) {
fprintf(stderr, "ERROR: failed to open file '%s'\n", opt->fn_chr_bin);
for (k = 0; k < kh_end(h); ++k)
if (kh_exist(h, k))
free((char*)kh_key(h, k));
cstr_ht_destroy(h);
return NULL;
}
MALLOC(ss, R_INF.total_reads); memset(ss, -1, sizeof((*ss))*R_INF.total_reads);
ks = ks_init(fp);
while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) {
khint_t k; ++n_tot;
phrase_hchar(str.s, &rname, &hid);
if((!(*rname)) || hid == (uint32_t)-1) {
fprintf(stderr, "ERROR: wrong hap id\n");
continue;
}
k = cstr_ht_get(h, rname);
if (k != kh_end(h)) {
ss[kh_val(h, k)] = hid;
if(hid > mhid) mhid = hid;
++n_bin;
// fprintf(stderr, "%s\t%u\trid::%ld\n", rname, hid, kh_val(h, k));
}
}
free(str.s);
ks_destroy(ks);
gzclose(fp);
for (k = 0; k < kh_end(h); ++k)
if (kh_exist(h, k))
free((char*)kh_key(h, k));
cstr_ht_destroy(h);
mhid++; CALLOC((*idx), mhid);
for (i = 0; i < (int64_t)R_INF.total_reads; ++i) {
if(ss[i] >= mhid) continue;
(*idx)[ss[i]] = 1;
}
*idx_n = mhid;
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> partitioned reads with external lists\n", __func__, yak_realtime(), yak_cpu_usage());
return ss;
}
void ha_triobin(const hifiasm_opt_t *opt)
{
memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads * sizeof(uint8_t));
anchor.cpp
+418 -28
View File
@@ -984,8 +984,8 @@ uint32_t *low_occ)
}
void minimizers_qgen0(ha_abuf_t *ab, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, Candidates_list *cl, kvec_t_u8_warp* k_flag,
void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ)
uint64_t minimizers_qgen0(ha_abuf_t *ab, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, Candidates_list *cl, kvec_t_u8_warp* k_flag,
void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint64_t ti_cut)
{
// fprintf(stderr, "+[M::%s]\n", __func__);
uint64_t i, k, l, max_cnt = UINT32_MAX, min_cnt = 0; int n, j; ha_mz1_t *z; seed1_t *s;
@@ -1042,7 +1042,7 @@ void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, kvec_t_u64_warp* dbg_ct, st_m
REALLOC(cl->list, cl->size);
}
k_mer_hit *p; uint64_t tid = (uint64_t)-1, tl = (uint64_t)-1;
k_mer_hit *p; uint64_t tid = (uint64_t)-1, tl = (uint64_t)-1, tcut_n = 0;
radix_sort_ha_an1(ab->a, ab->a + ab->n_a);
for (k = 1, l = 0; k <= ab->n_a; ++k) {
if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) {
@@ -1052,6 +1052,7 @@ void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, kvec_t_u64_warp* dbg_ct, st_m
tl = Get_READ_LENGTH((*rdb), tid);
// tl = rdb?Get_READ_LENGTH((*rdb), tid):udb->ug->u.a[tid].len;
}
if(tid < ti_cut) tcut_n = k;
for (i = l; i < k; i++) {
p = &cl->list[i];
p->readID = ab->a[i].srt>>33;
@@ -1078,6 +1079,7 @@ void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, kvec_t_u64_warp* dbg_ct, st_m
}
}
cl->length = ab->n_a;
return tcut_n;
}
void minimizers_qgen0_amz(ha_abuf_t *ab, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, Candidates_list *cl, kvec_t_u8_warp* k_flag,
@@ -1414,10 +1416,10 @@ uint32_t *low_occ, ha_mzl_t *in, uint64_t in_n, ha_mzl_t *idx, int64_t idx_n, ui
z = &in[i]; p = &(idx[z->x]); n = 0;
assert(z->pos == p->pos && z->rid == p->rid && z->span == p->span && z->rev == p->rev);
for (zi = z->x+1; zi < idx_n && idx[zi].x == p->x && n < mzl_cutoff; zi++) {
if(idx[zi].rid == rid) continue; n++;
if(idx[zi].rid == rid) {continue;} n++;
}
for (zi = ((int64_t)z->x)-1; zi >= 0 && idx[zi].x == p->x && n < mzl_cutoff; zi--) {
if(idx[zi].rid == rid) continue; n++;
if(idx[zi].rid == rid) {continue;} n++;
}
if((!n) || (n >= mzl_cutoff)) continue;
ab->n_a += n;
@@ -1435,11 +1437,11 @@ uint32_t *low_occ, ha_mzl_t *in, uint64_t in_n, ha_mzl_t *idx, int64_t idx_n, ui
///z is one of the minimizer
z = &in[i]; p = &(idx[z->x]); n = 0; ns = 0; ne = -1;
for (zi = z->x+1; zi < idx_n && idx[zi].x == p->x && n < mzl_cutoff; zi++) {
if(idx[zi].rid == rid) continue; n++;
if(idx[zi].rid == rid) {continue;} n++;
}
ne = zi;
for (zi = ((int64_t)z->x)-1; zi >= 0 && idx[zi].x == p->x && n < mzl_cutoff; zi--) {
if(idx[zi].rid == rid) continue; n++;
if(idx[zi].rid == rid) {continue;} n++;
}
ns = zi + 1;
if((!n) || (n >= mzl_cutoff)) continue;
@@ -1917,20 +1919,355 @@ void lchain_qgen_mcopy(Candidates_list* cl, overlap_region_alloc* ol, uint32_t r
for (i = 0; i < ol->length; ++i) ol->list[i].align_length = 0;
}
void lchain_qgen_mcopy_fast(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uint64_t rl, All_reads* rdb,
static inline uint64_t mz_pos_bsearch(const ha_mz1_t *a, uint64_t lo, uint64_t hi, uint64_t key) {
uint64_t mid;
if((lo >= hi) || (a[lo].pos > key)) {lo = 0;}
else if(a[lo].pos == key) {return lo;}
while (lo < hi) {
mid = lo + ((hi - lo) >> 1);
if (a[mid].pos < key) {
lo = mid + 1;
} else if(a[mid].pos > key) {
hi = mid;
} else {
lo = mid;
break;
}
}
if (lo < hi && a[lo].pos == key) {
return lo;
}
return ((uint64_t)-1);
}
uint8_t cmp_chain_aln(overlap_region *a, /**uint32_t ak0,**/ overlap_region *b, k_mer_hit *ca)
{
uint64_t ak, bk, al[2], bl[2]/**, af[2], bf[2]**/;
al[0] = a->non_homopolymer_errors; al[1] = a->overlapLen;
bl[0] = b->non_homopolymer_errors; bl[1] = b->overlapLen;
for (ak = al[0], bk = bl[0]; (ak < al[1]) && (bk < bl[1]) && (ca[ak].self_offset == ca[bk].self_offset); ak++, bk++);
// if(ka == 128 && kb == 129) {
// fprintf(stderr, "[M::%s::] al::[%lu,%lu), bl::[%lu,%lu), ak::%lu, bk::%lu\n", __func__, al[0], al[1], bl[0], bl[1], ak, bk);
// }
if((ak < al[1]) || (bk < bl[1])) return 0;
/**
af[0] = af[1] = bf[0] = bf[1] = ((uint64_t)-1);
for (ak = ak0, bk = bl[0]; (bk < bl[1]) && (ca[bk].self_offset < ca[ak].self_offset); bk++);
// if(!((bk < bl[1]) && (ca[ak].self_offset == ca[bk].self_offset))) {
// fprintf(stderr, "[M::%s::]\ta::%.*s->b::%.*s\n", __func__,
// (int32_t)Get_NAME_LENGTH(R_INF, a->y_id), Get_NAME(R_INF, a->y_id), (int32_t)Get_NAME_LENGTH(R_INF, b->y_id), Get_NAME(R_INF, b->y_id));
// }
assert((bk < bl[1]) && (ca[ak].self_offset == ca[bk].self_offset));
af[0] = ak; af[1] = ++ak;
bf[0] = bk; bf[1] = ++bk;
for (; (ak < al[1]) && (bk < bl[1]) && (ca[ak].self_offset == ca[bk].self_offset); ak++, bk++);
af[1] = ak; bf[1] = bk;
if(af[0] > al[0] && bf[0] > bl[0]) return 0;
if(af[1] < al[1] && bf[1] < bl[1]) return 0;
**/
return 1;
}
void debug_chain_aln_de(overlap_region_alloc *ol, Candidates_list *cl)
{
uint64_t k, z, zn; overlap_region *pk, *pz;
for (k = 0; k < ol->length; k++) {
pk = &(ol->list[k]);
for (z = zn = 0; z < ol->length; z++) {
if(z == k) continue;
pz = &(ol->list[z]);
if(cmp_chain_aln(pk, pz, cl->list)) {
if((pk->x_id != pz->x_id) || (pk->x_id == ((uint32_t)-1)) || (pz->x_id == ((uint32_t)-1))) {
fprintf(stderr, "[M::%s::type0] k::%lu, z::%lu\n", __func__, k, z);
exit(1);
}
zn++;
}
}
if((pk->x_id == ((uint32_t)-1)) && (zn > 0)) {
fprintf(stderr, "[M::%s::type1] k::%lu, zn::%lu\n", __func__, k, zn + 1);
exit(1);
}
}
}
void chain_aln_de(ha_abuf_t *ab, overlap_region_alloc *ol, Candidates_list *cl, asg32_v *ik)
{
uint64_t k, l, p, pi, c, z, zn, zk; uint32_t pn/**, ikn0 = ik->n**/;
// fprintf(stderr, "-0-[M::%s::] ik->n::%lu, ol->length::%lu\n", __func__, (uint64_t)ik->n, ol->length);
for (k = 0; k < ol->length; k++) {
// fprintf(stderr, "\n---[k::%lu::M::%s::%.*s(qid::%u)] q[%d, %d), t[%d, %d), sc::%d, occ::%u\n", k, __func__,
// (int32_t)Get_NAME_LENGTH(R_INF, ol->list[k].y_id), Get_NAME(R_INF, ol->list[k].y_id), ol->list[k].y_id,
// ol->list[k].x_pos_s, ol->list[k].x_pos_e+1, ol->list[k].y_pos_s, ol->list[k].y_pos_e+1, ol->list[k].shared_seed, ol->list[k].align_length);
if(ol->list[k].x_id != ((uint32_t)-1)) continue;
/**
for (l = ol->list[k].non_homopolymer_errors, zn = 0; l < ol->list[k].overlapLen; l++) {
p = cl->list[l].readID;
// pa = ik->a + (ab->mz.a[p].x>>32);
pi = (ab->mz.a[p].x>>32);
pn = ((uint32_t)(ab->mz.a[p].x));
assert(pn > 0);
fprintf(stderr, "[M::%s::] qpos::%u, tpos::%u, cnt::%u, pn::%u\n", __func__, cl->list[l].self_offset, cl->list[l].offset, cl->list[l].cnt>>8, pn);
for (z = 0; z < pn; z++) {
c = ik->a[pi + z];
// if(k == 128 && c == 129) fprintf(stderr, "+0+c::%lu\n", c);
if(c == k) continue;
// if(k == 128 && c == 129) fprintf(stderr, "+1+c::%lu\n", c);
if(ol->list[c].x_id != ((uint32_t)-1)) continue;
// if(k == 128 && c == 129) fprintf(stderr, "+2+c::%lu\n", c);
if(ol->list[c].strong) continue;
// if(k == 128 && c == 129) fprintf(stderr, "+3+c::%lu\n", c);
ol->list[c].strong = 1;
kv_push(uint32_t, *ik, c);
if(!cmp_chain_aln(k, &(ol->list[k]), l, c, &(ol->list[c]), cl->list)) continue;
// if(k == 128 && c == 129) fprintf(stderr, "+4+c::%lu\n", c);
ol->list[c].x_id = k;
zn++;
}
}
if(zn) ol->list[k].x_id = k;
for (l = ikn0; l < ik->n; l++) {
ol->list[ik->a[l]].strong = 0;
}
ik->n = ikn0;
**/
for (l = zk = ol->list[k].non_homopolymer_errors, zn = ((uint32_t)-1); l < ol->list[k].overlapLen; l++) {
p = cl->list[l].readID;
pn = ((uint32_t)(ab->mz.a[p].x));
assert(pn > 0);
if(pn < zn) {
zk = l; zn = pn;
}
// fprintf(stderr, "[M::%s::] qpos::%u, tpos::%u, cnt::%u, pn::%u\n", __func__, cl->list[l].self_offset, cl->list[l].offset, cl->list[l].cnt>>8, pn);
}
l = zk; p = cl->list[l].readID;
pi = (ab->mz.a[p].x>>32); pn = ((uint32_t)(ab->mz.a[p].x));
assert(pn > 0);
kv_push(uint32_t, *ik, k);
for (z = 0; z < pn; z++) {
c = ik->a[pi + z];
if(c <= k) continue;
if(ol->list[c].x_id != ((uint32_t)-1)) continue;
if(!cmp_chain_aln(&(ol->list[k]), &(ol->list[c]), cl->list)) continue;
ol->list[c].x_id = k;
zn++;
kv_push(uint32_t, *ik, c);
}
// if(zn) ol->list[k].x_id = k;
ol->list[k].x_id = k;
// fprintf(stderr, "[M::%s::] zn::%lu\n", __func__, zn);
}
/**
fprintf(stderr, "-1-[M::%s::] ik->n::%lu, ol->length::%lu\n", __func__, (uint64_t)ik->n, ol->length);
for (k = 0; k < ol->length; k++) {
if((ol->list[k].x_id != ((uint32_t)-1)) && (ol->list[k].x_id != k)) continue;
if(k == ol->list[k].x_id) {
fprintf(stderr, "[M::%s::] cs::%lu->", __func__, k);
for (l = zn = 0; l < ol->length; l++) {
if(ol->list[k].x_id == ol->list[l].x_id) {
fprintf(stderr, "%lu,", l);
zn++;
}
}
fprintf(stderr, "(znn::%lu)\n", zn);
} else {
fprintf(stderr, "[M::%s::] cs::%lu->(znn::1)\n", __func__, k);
}
for (l = ol->list[k].non_homopolymer_errors, zn = 0; l < ol->list[k].overlapLen; l++) {
fprintf(stderr, "[M::%s::] qpos::%u, tpos::%u, cnt::%u\n", __func__, cl->list[l].self_offset, cl->list[l].offset, cl->list[l].cnt>>8);
}
}
debug_chain_aln_de(ol, cl);
// for (k = ik->n - ol->length; k < ik->n; k++) {
// ik->a[k]
// }
**/
}
void gen_chain_clus(ha_abuf_t *ab, overlap_region_alloc *ol, Candidates_list *cl, asg32_v *ik)
{
if(ol->length <= 0) return;
ks_introsort_or_ss(ol->length, ol->list);
uint64_t k, l, c, p_b, p_r, cln = cl->length; uint32_t *pa = NULL, pn, xid0 = ol->list[0].x_id;
for (k = 1, l = 0; k < ol->length; k++) {///sort for merging
if ((k == ol->length) || (ol->list[l].shared_seed != ol->list[k].shared_seed)) {
if (k - l > 1) {
ks_introsort_or_xs(k - l, ol->list + l);
}
l = k;
}
}
for (k = 0; k < ab->mz.n; k++) ab->mz.a[k].x = 0;
for (k = 0; k < ol->length; k++) {
l = ol->list[k].non_homopolymer_errors;
p_r = cl->list[l].readID;
for (; (l < cln) && (cl->list[l].readID == p_r); l++) {;}
ol->list[k].overlapLen = l;
}
for (k = ik->n = 0, p_b = ab->mz.n; k < ol->length; k++) {
for (l = ol->list[k].non_homopolymer_errors; l < ol->list[k].overlapLen; l++) {
p_b = mz_pos_bsearch(ab->mz.a, p_b + 1, ab->mz.n, cl->list[l].self_offset);
assert(p_b != ((uint64_t)-1));
// idx->a[p]++;
ab->mz.a[p_b].x++; ik->n++;
cl->list[l].readID = p_b;///minimizer id
}
// fprintf(stderr, "[M_beg::%s::]\tk::%lu\ty_id::%u\tc_n::%u\tc_beg::%u\n", __func__, k, ol->list[k].y_id, ol->list[k].overlapLen - ol->list[k].non_homopolymer_errors, ol->list[k].non_homopolymer_errors);
}
/**
for (k = ik->n = 0, p_b = ab->mz.n; k < ol->length; k++) {
l = ol->list[k].non_homopolymer_errors;
p_r = cl->list[l].readID;
for (; (l < cln) && (cl->list[l].readID == p_r); l++) {
p_b = mz_pos_bsearch(ab->mz.a, p_b + 1, ab->mz.n, cl->list[l].self_offset);
assert(p_b != ((uint64_t)-1));
// idx->a[p]++;
ab->mz.a[p_b].x++; ik->n++;
cl->list[l].readID = p_b;///minimizer id
}
ol->list[k].overlapLen = l;
fprintf(stderr, "[M_beg::%s::]\tk::%lu\ty_id::%u\tc_n::%u\tc_beg::%u\n", __func__, k, ol->list[k].y_id, ol->list[k].overlapLen - ol->list[k].non_homopolymer_errors, ol->list[k].non_homopolymer_errors);
}
**/
kv_resize(uint32_t, *ik, ik->n);
memset(ik->a, -1, sizeof((*(ik->a)))*ik->n);
for (k = l = 0; k < ab->mz.n; k++) {
ab->mz.a[k].x |= (l<<32);
l += (uint32_t)(ab->mz.a[k].x);
if((uint32_t)(ab->mz.a[k].x)) {
ik->a[(ab->mz.a[k].x>>32) + ((uint32_t)(ab->mz.a[k].x)) - 1] = 0;
}
}
assert(l == ik->n);
for (k = 0; k < ol->length; k++) {
for (l = ol->list[k].non_homopolymer_errors; l < ol->list[k].overlapLen; l++) {
p_b = cl->list[l].readID;
pa = ik->a + (ab->mz.a[p_b].x>>32);
pn = ((uint32_t)(ab->mz.a[p_b].x));
assert(pn > 0);
c = pa[pn - 1];
assert(c < pn);
pa[c] = k; ///pa[c] = l;
if(c + 1 < pn) {pa[pn - 1]++;}
// cl->list[l].readID = k;
}
// ol->list[k].overlapLen = 0;
ol->list[k].x_id = ((uint32_t)-1);///mark for clustering
}
chain_aln_de(ab, ol, cl, ik);
for (k = ik->n - ol->length, l = 0; k < ik->n; k++) {
ik->a[l++] = ik->a[k];
}
ik->n = ol->length;
ik->n = ol->length<<1; kv_resize(uint32_t, *ik, ik->n);
memset(ik->a + ol->length, -1, (sizeof((*(ik->a)))*ol->length));
// ik->n = ol->length;
for (k = 1, l = 0; k <= ol->length; k++) {
if((k == ol->length) || (ol->list[ik->a[l]].x_id != ol->list[ik->a[k]].x_id)) {
for (c = l; c < k; c++) {
ik->a[ik->a[c]+ol->length] = l;
}
l = k;
}
}
// for (k = 0, ik->n = ol->length; k < ol->length; k++) {
// ik->a[ik->n++] = ((ol->list[k].x_id == ((uint32_t)-1))?k:ol->list[k].x_id);
// }
///reset
for (k = 0; k < ol->length; k++) {
// fprintf(stderr, "[M_end::%s::]\tk::%lu\ty_id::%u\tc_n::%u\tc_beg::%u\n", __func__, k, ol->list[k].y_id, ol->list[k].overlapLen - ol->list[k].non_homopolymer_errors, ol->list[k].non_homopolymer_errors);
for (l = ol->list[k].non_homopolymer_errors; l < ol->list[k].overlapLen; l++) {
cl->list[l].readID = k;
}
ol->list[k].overlapLen = 0;
ol->list[k].x_id = xid0;
}
}
void lchain_qgen_mcopy_fast(ha_abuf_t *ab, Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uint64_t rl, All_reads* rdb,
uint32_t apend_be, uint64_t max_n_chain, 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 quick_check,
uint32_t gen_off, int64_t mcopy_num, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, st_mt_t *sp, uint64_t ocv_w)
int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate_h, uint64_t cl_hn, double bw_rate_l, uint64_t cl_ln, int64_t quick_check,
uint32_t gen_off, int64_t mcopy_num, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, st_mt_t *sp, uint64_t ocv_w, uint8_t is_raw_chain)
{
// fprintf(stderr, "+[M::%s] chain_cutoff::%u\n", __func__, chain_cutoff);
uint64_t i, k, l, m, cn = cl->length, yid, ol0, lch, *cc = NULL, cwn = 0, cws, cwe, os, oe, rs, re, cw0, cw1/**, dbgn = 0**/; overlap_region *r, t; ///srt = 0
clear_overlap_region_alloc(ol);
// for (l = 0, k = 1, m = 0, lch = 0; k <= cn; k++) {
// if((k == cn) || (cl->list[k].readID != cl->list[l].readID)) {
// if(cl->list[l].readID != rid) {
// yid = cl->list[l].readID; ol0 = ol->length;
// m += lchain_qdp_mcopy_fast(cl, l, k-l, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate,
// rid, rl, Get_READ_LENGTH((*rdb), yid), quick_check, apend_be, gen_off, mcopy_num, mcopy_rate, mcopy_khit_cut, 1);
// if((chain_cutoff >= 2) && (!lch)) {
// for (i = ol0; (i<ol->length) && (!lch); i++) {
// if(ol->list[i].align_length < chain_cutoff) lch = 1;
// }
// }
// }
// l = k;
// }
// }
// cl->length = m;
///cl->list[0, cl_hn) & cl->list[cl_hn, cl_hn + cl_ln)
cn = cl_hn;
for (l = 0, k = 1, m = 0, lch = 0; k <= cn; k++) {
if((k == cn) || (cl->list[k].readID != cl->list[l].readID)) {
if(cl->list[l].readID != rid) {
yid = cl->list[l].readID; ol0 = ol->length;
m += lchain_qdp_mcopy_fast(cl, l, k-l, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate,
m += lchain_qdp_mcopy_fast(cl, l, k-l, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate_h,
rid, rl, Get_READ_LENGTH((*rdb), yid), quick_check, apend_be, gen_off, mcopy_num, mcopy_rate, mcopy_khit_cut, 1);
if((chain_cutoff >= 2) && (!lch)) {
for (i = ol0; (i<ol->length) && (!lch); i++) {
@@ -1941,14 +2278,47 @@ void lchain_qgen_mcopy_fast(Candidates_list* cl, overlap_region_alloc* ol, uint3
l = k;
}
}
cn = cl_hn + cl_ln;
for (; k <= cn; k++) {
if((k == cn) || (cl->list[k].readID != cl->list[l].readID)) {
if(cl->list[l].readID != rid) {
yid = cl->list[l].readID; ol0 = ol->length;
m += lchain_qdp_mcopy_fast(cl, l, k-l, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate_l,
rid, rl, Get_READ_LENGTH((*rdb), yid), quick_check, apend_be, gen_off, mcopy_num, mcopy_rate, mcopy_khit_cut, 1);
if((chain_cutoff >= 2) && (!lch)) {
for (i = ol0; (i<ol->length) && (!lch); i++) {
if(ol->list[i].align_length < chain_cutoff) lch = 1;
}
}
}
l = k;
}
}
cl->length = m;
// fprintf(stderr, "[M::%s::] rn::%lu\tmax_n_chain::%lu\n", __func__, ol->length, max_n_chain);
// for (k = 0; k < ol->length; k++) {
// fprintf(stderr, "---[M::%s::%.*s(qid::%u)] q[%d, %d), t[%d, %d), sc::%d, type::%d\n", __func__,
// (int32_t)Get_NAME_LENGTH(R_INF, ol->list[k].y_id), Get_NAME(R_INF, ol->list[k].y_id), ol->list[k].y_id,
// ol->list[k].x_pos_s, ol->list[k].x_pos_e+1, ol->list[k].y_pos_s, ol->list[k].y_pos_e+1, ol->list[k].shared_seed, ha_ov_type(&(ol->list[k]), rl));
// }
// fprintf(stderr, "\n[M::%s::] rn::%lu\tmax_n_chain::%lu\n", __func__, ol->length, max_n_chain);
/**
for (k = 0; k < ol->length; k++) {
if(ha_ov_type(&(ol->list[k]), rl) == 1) {
fprintf(stderr, "---[M::%s::%.*s(qid::%u)] q[%d, %d), t[%d, %d), sc::%d, occ::%u, type::%d\n", __func__,
(int32_t)Get_NAME_LENGTH(R_INF, ol->list[k].y_id), Get_NAME(R_INF, ol->list[k].y_id), ol->list[k].y_id,
ol->list[k].x_pos_s, ol->list[k].x_pos_e+1, ol->list[k].y_pos_s, ol->list[k].y_pos_e+1, ol->list[k].shared_seed, ol->list[k].align_length, ha_ov_type(&(ol->list[k]), rl));
int64_t km;
for (km = ol->list[k].non_homopolymer_errors; (km < cl->length) && (cl->list[km].readID == cl->list[ol->list[k].non_homopolymer_errors].readID); km++) {
fprintf(stderr, "[M::%s::] qpos::%u, tpos::%u, cnt::%u\n", __func__, cl->list[km].self_offset, cl->list[km].offset, cl->list[km].cnt>>8);
}
fprintf(stderr, "[M::%s::]\tk::%lu\tsi::%u\tei::%ld\n", __func__, k, ol->list[k].non_homopolymer_errors, km);
}
}
**/
// gen_chain_clus(ab, ol, cl, v32);
if(is_raw_chain) return;
k = ol->length;
if (ol->length > max_n_chain) {
@@ -1961,7 +2331,7 @@ void lchain_qgen_mcopy_fast(Candidates_list* cl, overlap_region_alloc* ol, uint3
++n[w];
if (((uint64_t)n[w]) == max_n_chain) s[w] = r->shared_seed;
}
// fprintf(stderr, "[M::%s::] s[0]::%d, s[1]::%d, s[2]::%d, s[3]::%d\n", __func__, s[0], s[1], s[2], s[3]);
/**fprintf(stderr, "top[M::%s::] s[0]::%d, s[1]::%d, s[2]::%d, s[3]::%d\n", __func__, s[0], s[1], s[2], s[3]);**/
if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
if((((uint64_t)n[3]) >= max_n_chain) && (rl >= ocv_w)) {
cwn = (rl/ocv_w) + ((rl % ocv_w)?(1):(0));
@@ -1972,7 +2342,7 @@ void lchain_qgen_mcopy_fast(Candidates_list* cl, overlap_region_alloc* ol, uint3
assert(cwe > cws);
cc[i] = (cwe - cws)*(max_n_chain>>1);
// fprintf(stderr, "[M::%s::] cw::[%lu, %lu), cc::%lu\n", __func__, cws, cwe, cc[i]);
if(cc[i] > UINT32_MAX) cc[i] = UINT32_MAX; cc[i] <<= 32;
if(cc[i] > UINT32_MAX) {cc[i] = UINT32_MAX;} cc[i] <<= 32;
cws += ocv_w;
}
}
@@ -2099,6 +2469,11 @@ void lchain_qgen_mcopy_fast(Candidates_list* cl, overlap_region_alloc* ol, uint3
// fprintf(stderr, "+[M::%s] rid::%u, ol->length0::%lu, dbgn::%lu\n", __func__, rid, ol->length, dbgn);
}
void srt_olst(overlap_region_alloc* ol)
{
ks_introsort_or_xs(ol->length, ol->list);
}
void lchain_qgen_mcopy_fast_re1(Candidates_list* cl, uint32_t cl_beg, overlap_region_alloc* ol, uint32_t rid, uint64_t rl, uint64_t tl,
uint32_t apend_be, 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 quick_check,
@@ -2300,18 +2675,33 @@ void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t
}
void h_ec_lchain(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, All_reads *rref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres,
int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, 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, int64_t mcopy_num, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, uint64_t ocv_w)
int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, 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, int64_t mcopy_num, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, uint64_t ocv_w, uint8_t is_raw_chain)
{
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
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);
// minimizers_gen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, dbg_ct, sp, high_occ, low_occ);
minimizers_qgen0(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, rref, dbg_ct, sp, high_occ, low_occ);
minimizers_qgen0(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, rref, dbg_ct, sp, high_occ, low_occ, ((uint64_t)-1));
// lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
// lchain_qgen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
///no need to sort here, overlap_list has been sorted at lchain_gen
lchain_qgen_mcopy_fast(cl, overlap_list, rid, rl, rref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off, mcopy_num, mcopy_rate, chain_cutoff, mcopy_khit_cut, sp, ocv_w);
lchain_qgen_mcopy_fast(ab, cl, overlap_list, rid, rl, rref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, cl->length, bw_thres, 0, quick_check, gen_off, mcopy_num, mcopy_rate, chain_cutoff, mcopy_khit_cut, sp, ocv_w, is_raw_chain);
}
void h_ec_lchain_hybrid(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, All_reads *rref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres_h, double bw_thres_l,
int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, 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, int64_t mcopy_num, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, uint64_t ocv_w, uint64_t ti_cut, uint8_t is_raw_chain)
{
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip; uint64_t tcut_n = 0;
set_lchain_dp_op(is_accurate, mz_k, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check);
// minimizers_gen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, dbg_ct, sp, high_occ, low_occ);
tcut_n = minimizers_qgen0(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, rref, dbg_ct, sp, high_occ, low_occ, ti_cut);
// lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
// lchain_qgen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
///no need to sort here, overlap_list has been sorted at lchain_gen
lchain_qgen_mcopy_fast(ab, cl, overlap_list, rid, rl, rref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres_h, tcut_n, bw_thres_l, cl->length - tcut_n, quick_check, gen_off, mcopy_num, mcopy_rate, chain_cutoff, mcopy_khit_cut, sp, ocv_w, is_raw_chain);
}
void h_ec_lchain_amz(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, All_reads *rref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres,
@@ -2326,7 +2716,7 @@ void h_ec_lchain_amz(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_
// lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
// lchain_qgen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
///no need to sort here, overlap_list has been sorted at lchain_gen
lchain_qgen_mcopy_fast(cl, overlap_list, rid, rl, rref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off, enable_mcopy, mcopy_rate, chain_cutoff, mcopy_khit_cut, sp, ocv_w);
lchain_qgen_mcopy_fast(ab, cl, overlap_list, rid, rl, rref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, cl->length, bw_thres, 0, quick_check, gen_off, enable_mcopy, mcopy_rate, chain_cutoff, mcopy_khit_cut, sp, ocv_w, 0);
}
uint64_t recalu_minimizer0(char *s, uint64_t len, uint64_t is_hpc, int64_t mz_k, uint64_t mz_h, tiny_queue_t *tq, uint64_t *rpos, uint64_t *rspan)
@@ -3050,10 +3440,10 @@ uint64_t recalu_minimizer_non_retrieve(uint64_t rid, anchor1_t *z, asg16_v *sc,
}
c = tstr[s1];
for (k = s1 - 1; k >= 0 && tstr[k] == c; k--); s1 = k + 1;
for (k = s1 - 1; k >= 0 && tstr[k] == c; k--){;} s1 = k + 1;
c = tstr[e1-1];
for (k = e1; k < tl && tstr[k] == c; k++); e1 = k;
for (k = e1; k < tl && tstr[k] == c; k++){;} e1 = k;
if(!trev) {
s0 = s1; e0 = e1;
@@ -3075,15 +3465,15 @@ uint64_t recalu_minimizer_non_retrieve(uint64_t rid, anchor1_t *z, asg16_v *sc,
void hpc_ext_check(All_reads *rref, uint32_t id, int64_t s0, int64_t e0, int64_t l, uint64_t rev, char *buf)
{
if(s0 < 0) s0 = 0; if(e0 > l) e0 = l;
if(s0 < 0) {s0 = 0;} if(e0 > l) {e0 = l;}
int64_t s = s0 - 256, e = e0 + 256, n, k, os0, oe0, os1, oe1; char c; if(s < 0) s = 0; if(e > l) e = l;
recover_UC_Read_sub_region(buf, s, e - s, rev, rref, id);
os0 = s0 - s; oe0 = e0 - s; n = e - s; os1 = os0; oe1 = oe0;
c = buf[os0];
for (k = os0 - 1; k >= 0 && buf[k] == c; k--); os1 = k + 1;
for (k = os0 - 1; k >= 0 && buf[k] == c; k--){;} os1 = k + 1;
c = buf[oe0-1];
for (k = oe0; k < n && buf[k] == c; k++); oe1 = k;
for (k = oe0; k < n && buf[k] == c; k++){;} oe1 = k;
if(os1 != os0 || oe1 != oe0) {
fprintf(stderr, "[M::%s]\to0::[%ld,%ld)\to1::[%ld,%ld)\n", __func__, os0 + s0, oe0 + s0, os1 + s0, oe1 + s0);
ecovlp.cpp
+846 -155
View File
File diff suppressed because it is too large Load Diff
gfa_ut.cpp
+69 -51
View File
@@ -966,7 +966,7 @@ void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max
{
asg64_v tx = {0,0,0}, *b = NULL;
uint32_t v, w, i, k, n_vtx = g->n_seq<<1;
asg_arc_t *av, *aw, *ve, *vmax, *we; uint32_t nv, nw, kv, kw, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou;
asg_arc_t *av = NULL, *aw = NULL, *ve = NULL, *vmax, *we = NULL; uint32_t nv, nw, kv, kw, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou;
uint32_t trioF = (uint32_t)-1, ntrioF = (uint32_t)-1;
if(in) b = in;
else b = &tx;
@@ -1280,7 +1280,7 @@ uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, ma_hit_t_alloc *rev, R_to_
{
asg64_v tx = {0,0,0}, *b = NULL;
uint32_t i, k, v, w, n_vtx = g->n_seq<<1, nv, nw, kv, kw, trioF = (uint32_t)-1, ntrioF = (uint32_t)-1, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou;
asg_arc_t *av, *aw, *ve, *we, *vl_max, *wl_max;
asg_arc_t *av = NULL, *aw = NULL, *ve = NULL, *we = NULL, *vl_max = NULL, *wl_max = NULL;
if(in) b = in;
else b = &tx;
@@ -1429,7 +1429,7 @@ uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, uint32_t test_bub, ma_hit_
{
asg64_v tx = {0,0,0}, *b = NULL;
uint32_t i, k, v, w, wz, n_vtx = g->n_seq<<1, nv, nw, kv, kw, trioF = (uint32_t)-1, ntrioF = (uint32_t)-1, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou, olw[2], m;
asg_arc_t *av, *aw, *ve, *we, *vl_max, *wl_max; ma_hit_t_alloc *z = NULL;
asg_arc_t *av = NULL, *aw = NULL, *ve = NULL, *we = NULL, *vl_max = NULL, *wl_max = NULL; ma_hit_t_alloc *z = NULL;
if(in) b = in;
else b = &tx;
@@ -1576,7 +1576,7 @@ uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, uint32_t test_bub, ma_hit_
void asg_arc_cut_chimeric_cmk(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_topo, uint32_t min_ou, uint32_t test_bub, uint8_t *cmk, uint32_t cmk_cut)
{
asg64_v tx = {0,0,0}, *b = NULL; asg_arc_t *av, *aw, *ve, *we;
asg64_v tx = {0,0,0}, *b = NULL; asg_arc_t *av = NULL, *aw = NULL, *ve = NULL, *we = NULL;
uint32_t i, k, v, w, n_vtx = g->n_seq<<1, nv, nw, kv, kw, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou;
if(in) b = in;
@@ -1681,7 +1681,7 @@ uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *
{
asg64_v tx = {0,0,0}, *b = NULL;
uint32_t i, k, v, w, n_vtx = g->n_seq<<1, nv, nw, kv, kw, trioF = (uint32_t)-1, ntrioF = (uint32_t)-1, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou;
asg_arc_t *av, *aw, *ve, *we, *vl_max, *wl_max;
asg_arc_t *av = NULL, *aw = NULL, *ve = NULL, *we = NULL, *vl_max = NULL, *wl_max = NULL;
if(in) b = in;
else b = &tx;
@@ -2247,7 +2247,7 @@ void asg_arc_cut_contain(flex_asg_t *fg, asg64_v *in, asg64_v *in0, R_to_U* rI,
// }
}
// stats_sysm(g);
if(!in) free(tx.a); if(!in0) free(tx0.a);
if(!in) {free(tx.a);} if(!in0) {free(tx0.a);}
if(cnt > 0) flex_asg_t_cleanup(fg);
// fprintf(stderr, "-[M::%s]\n", __func__);
}
@@ -2936,9 +2936,9 @@ void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd)
for (k = 0; k < nv; ++k) {
if ((av[k].v>>1) == dst) {
w = av[k].v;
fprintf(stderr, "[M::%s::]\t%.*s(%c)\t%.*s(%c)\tou::%u\n", __func__,
fprintf(stderr, "[M::%s::]\t%.*s(%c)\t%.*s(%c)\tou::%u\tdel::%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[k].ou);
(int)Get_NAME_LENGTH(R_INF, (w>>1)), Get_NAME(R_INF, (w>>1)), "+-"[w&1], av[k].ou, av[k].del);
}
}
@@ -2947,9 +2947,9 @@ void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd)
for (k = 0; k < nv; ++k) {
if ((av[k].v>>1) == dst) {
w = av[k].v;
fprintf(stderr, "[M::%s::]\t%.*s(%c)\t%.*s(%c)\tou::%u\n", __func__,
fprintf(stderr, "[M::%s::]\t%.*s(%c)\t%.*s(%c)\tou::%u\tdel::%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[k].ou);
(int)Get_NAME_LENGTH(R_INF, (w>>1)), Get_NAME(R_INF, (w>>1)), "+-"[w&1], av[k].ou, av[k].del);
}
}
}
@@ -3036,6 +3036,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
double drop = min_ovlp_drop_ratio;
int64_t i; asg64_v bu = {0,0,0}, ba = {0,0,0}; uint32_t l_drop = 2000; flex_asg_t *fg = NULL; uint32_t min_diff = 0, step_diff = 2000;
if(is_ou) fg = init_flex_asg_t(sg, uopt->sources, uopt->min_ovlp, uopt->max_hang, asm_opt.max_hang_rate, gap_fuzz);
// prt_specfic_sge(sg, 22708, 22646, "--sa--");
// if(is_ou) update_sg_uo(sg, src);///do not do it here
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
// exit(1);
@@ -3071,7 +3072,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
}
// fprintf(stderr, "(0):i->%ld, drop->%f\n", i, drop);
// prt_specfic_sge(sg, 10531, 10519, "--0--");
// prt_specfic_sge(sg, 22708, 22646, "--0--");
// print_vw_edge(sg, 34156, 34090, "0");
// stats_chimeric(sg, src, &bu);
@@ -3080,18 +3081,23 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
asg_arc_cut_chimeric(sg, src, &bu, is_ou?ou_thres:(uint32_t)-1, uopt->te);///p_telo
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// prt_specfic_sge(sg, 10531, 10519, "--1--");
// 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**/);
// debug_edges(&dbg, d, 2);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// prt_specfic_sge(sg, 10531, 10519, "--2--");
// 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);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// prt_specfic_sge(sg, 10531, 10519, "--3--");
// 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);
// // exit(1);
// }
// prt_specfic_sge(sg, 22708, 22646, "--3--");
// if(is_ou) asg_arc_cut_contain(fg, &bu, &ba, rI, ((i+1)<clean_round)?ou_drop_rate:-1);
if(is_ou) {
asg_arc_cut_contain(fg, &bu, &ba, rI, ou_drop_rate, 0);
@@ -3100,17 +3106,13 @@ 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);
// prt_specfic_sge(sg, 10531, 10519, "--4--");
// 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);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// prt_specfic_sge(sg, 10531, 10519, "--5--");
// prt_specfic_sge(sg, 22708, 22646, "--5--");
// if(i == 3) {
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty3.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
// // exit(1);
// }
/**
if(cmk && asm_opt.chemical_cov > FORCE_CUT) {
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0);
@@ -3126,6 +3128,8 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
}
}
// prt_specfic_sge(sg, 22708, 22646, "--sb--");
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);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
@@ -3133,7 +3137,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
if(is_ou) min_diff = step_diff;
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty4.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
// debug_info_of_specfic_node("m64012_190921_234837/111673711/ccs", sg, rI, "end");
// debug_info_of_specfic_node("7b70a587-f56c-48ac-adf8-b98a67365063_2", sg, rI, "end");
// debug_info_of_specfic_node("m64011_190830_220126/95028102/ccs", sg, rI, "end");
if(is_ou) {
asg_arc_cut_contain(fg, &bu, &ba, rI, ou_drop_rate, 0);
@@ -3143,12 +3147,17 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
}
if(!is_ou) asg_iterative_semi_circ(sg, src, &bu, max_tip, 1, uopt->te);
// prt_specfic_sge(sg, 22708, 22646, "--sb-0---");
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty5.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_cut_large_indel(sg, &bu, max_tip, HARD_OU_DROP, is_ou, min_diff);///shoule we ignore ou here?
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
// 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);
if(!is_ou) {
@@ -3169,9 +3178,13 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
}
// prt_specfic_sge(sg, 22708, 22646, "--sb-2---");
// 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);
// prt_specfic_sge(sg, 22708, 22646, "--sb-3---");
/**
rescue_contained_reads_aggressive(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 10, 1, 0, NULL, NULL, b_mask_t);
rescue_missing_overlaps_aggressive(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 1, 0, NULL, b_mask_t);
@@ -3185,9 +3198,13 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
**/
post_rescue(uopt, sg, src, rev, rI, b_mask_t, is_ou, cmk);
// prt_specfic_sge(sg, 22708, 22646, "--sb-4---");
ug_ext_gfa(uopt, sg, ug_ext_len);
// if(is_ou) dedup_contain_g(uopt, sg);
// prt_specfic_sge(sg, 22708, 22646, "--sb-5---");
// exit(1)
output_unitig_graph(sg, uopt->coverage_cut, o_file, src, rI, uopt->max_hang, uopt->min_ovlp);
@@ -3203,6 +3220,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
if(is_ou) {
des_flex_asg_t(fg); free(fg);
}
// prt_specfic_sge(sg, 22708, 22646, "--sb-6---");
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
// exit(1);
// print_node(sg, 17078); //print_node(sg, 8311); print_node(sg, 8294);
@@ -6065,7 +6083,7 @@ void poa_chain_0(poa_g_t *g, ma_ug_t *ug, uc_block_t *raw, ul_str_t *str, int64_
i = integer_g_chain(g, ug, buf->b.a, buf->b.n, buf, &rr); assert(i);
if(!i) return; buf->b.n = rr.e; assert(buf->b.n);
if(!i) {return;} buf->b.n = rr.e; assert(buf->b.n);
append_aligned_integer_seq_by_aln_pair(ug, raw, g, g->seq.n, pat, pat_n, is_rev, buf->b.a, buf->b.n, str_id, (is_rev?(str->cn - e):(s)));
(*is_circle) = update_poa_dp(g, debug_qid, ug);
}
@@ -8741,7 +8759,7 @@ uint32_t ulg_arc_cut_tips(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t max_ext, uin
}
// stats_sysm(g);
if(!in) free(tx.a); if(!ib) free(tb.a);
if(!in) {free(tx.a);} if(!ib) {free(tb.a);}
if (cnt > 0) asg_cleanup(g);
return cnt;
@@ -8817,7 +8835,7 @@ float len_rat, uint32_t is_trio, uint32_t topo_level, uint32_t hom_check, uint32
{
asg64_v tx = {0,0,0}, tb = {0,0,0}, *b = NULL, *ub = NULL; asg_t *g = ug->g;
uint32_t v, w, i, k, kv, kw, nv, nw, cnt = 0, n_vtx = g->n_seq<<1, n_het, n_hom, ff, ol_max, mm_ol, to_del;
asg_arc_t *av, *aw, *ve, *we, *vl_max, *wl_max;
asg_arc_t *av = NULL, *aw = NULL, *ve = NULL, *we = NULL, *vl_max = NULL, *wl_max = NULL;
b = (in?(in):(&tx)); ub = (ib?(ib):(&tb));
for (v = b->n = 0; v < n_vtx; ++v) {
@@ -8923,7 +8941,7 @@ float len_rat, uint32_t is_trio, uint32_t topo_level, uint32_t hom_check, uint32
}
if(!in) free(tx.a); if(!ib) free(tb.a);
if(!in) {free(tx.a);} if(!ib) {free(tb.a);}
if (cnt > 0) asg_cleanup(g);
return cnt;
@@ -8994,7 +9012,7 @@ uint32_t is_trio, uint32_t topo_level, asg64_v *in, asg64_v *ib)
}
}
if(!in) free(tx.a); if(!ib) free(tb.a);
if(!in) {free(tx.a);} if(!ib) {free(tb.a);}
if (cnt > 0) asg_cleanup(g);
return cnt;
@@ -9006,7 +9024,7 @@ uint32_t *ul_cnt, uint32_t *bridge_w, uint32_t *bridge_am, asg64_v *b)
{
uint32_t v = s, w = (uint32_t)-1, kv, kw, uv = (uint32_t)-1, uw = (uint32_t)-1, bv, bw, l;
ul_bg_t *bg = &(uidx->uovl.bg);
if(occ) (*occ) = 0; if(ul_cnt) (*ul_cnt) = 0; if(bridge_w) (*bridge_w) = 0; if(bridge_am) (*bridge_am) = 0;
if(occ) {(*occ) = 0;} if(ul_cnt) {(*ul_cnt) = 0;} if(bridge_w) {(*bridge_w) = 0;} if(bridge_am) {(*bridge_am) = 0;}
while (1) {
if(occ) (*occ)++;
@@ -9181,7 +9199,7 @@ float len_rat, uint32_t is_trio, uint32_t topo_level, uint32_t *max_drop_len, as
}
if(!in) free(tx.a); if(!ib) free(tb.a);
if(!in) {free(tx.a);} if(!ib) {free(tb.a);}
if (cnt > 0) asg_cleanup(g);
return cnt;
}
@@ -9453,7 +9471,7 @@ float len_rat, uint32_t is_trio, uint32_t topo_level, uint32_t *max_drop_len, as
if (nv < 2) continue;
for (i = kv = 0, ul_max = -1; i < nv; ++i) {
if(av[i].del) continue; kv++;
if(av[i].del) {continue;} kv++;
get_ul_path_info(uidx, ug, av[i].v, NULL, NULL, &ul_cnt, NULL, NULL, NULL);
if(ul_cnt > ul_max) ul_max = ul_cnt;
}
@@ -9546,7 +9564,7 @@ float len_rat, uint32_t is_trio, uint32_t topo_level, uint32_t *max_drop_len, as
}
if(!in) free(tx.a); if(!ib) free(tb.a);
if(!in) {free(tx.a);} if(!ib) {free(tb.a);}
if (cnt > 0) asg_cleanup(g);
return cnt;
}
@@ -9774,7 +9792,7 @@ asg64_v *in, asg64_v *ib)
if (nv < 2) continue;
for (i = kv = 0; i < nv && kv < 2; ++i) {
if(av[i].del) continue; kv++;
if(av[i].del) {continue;} kv++;
}
if(kv < 2) continue;
@@ -9877,7 +9895,7 @@ asg64_v *in, asg64_v *ib)
}
}
if(!in) free(tx.a); if(!ib) free(tb.a);
if(!in) {free(tx.a);} if(!ib) {free(tb.a);}
if (cnt > 0) asg_cleanup(g);
return cnt;
}
@@ -10134,7 +10152,7 @@ uint64_t ulg_pop_bubble(ul_resolve_t *uidx, ma_ug_t *ug, uint64_t* i_max_dist, u
free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
if(n_pop) asg_cleanup(g);
if(!in) free(tx.a); if(!ib) free(tb.a);
if(!in) {free(tx.a);} if(!ib) {free(tb.a);}
// fprintf(stderr, "[M::%s::] Done...\n", __func__);
return n_pop;
}
@@ -10348,7 +10366,7 @@ uint32_t skip_hom, uint32_t *max_drop_len, asg64_v *in, asg64_v *ib)
av = asg_arc_a(g, v); nv = asg_arc_n(g, v);
if (nv < 2) continue;
for (i = kv = 0; i < nv && kv <= 2; ++i) {
if(av[i].del) continue; kv++;
if(av[i].del) {continue;} kv++;
}
if(kv != 2) continue;
for (i = 0; i < nv; ++i) {
@@ -10497,7 +10515,7 @@ uint32_t skip_hom, uint32_t *max_drop_len, asg64_v *in, asg64_v *ib)
if(!in) free(tx.a); if(!ib) free(tb.a);
if(!in) {free(tx.a);} if(!ib) {free(tb.a);}
if (cnt > 0) asg_cleanup(g);
// fprintf(stderr, "[M::%s::] cnt::%u\n", __func__, cnt);
return cnt;
@@ -10971,7 +10989,7 @@ uint32_t is_topo, uint32_t *max_drop_len)
// }
asg64_v tx = {0,0,0}, tz = {0,0,0}, *b = NULL, *ub = NULL;
uint32_t i, k, v, w, n_vtx = g->n<<1, nv, nw, kv, kw, /**trioF = (uint32_t)-1, ntrioF = (uint32_t)-1,**/ ol_max, ou_max, to_del, cnt = 0, mm_ol;
usg_arc_t *av, *aw, *ve, *we; uint64_t x, kocc[2], ou; uint8_t *f; CALLOC(f, g->n);
usg_arc_t *av = NULL, *aw = NULL, *ve = NULL, *we = NULL; uint64_t x, kocc[2], ou; uint8_t *f; CALLOC(f, g->n);
b = ((in_0)?(in_0):(&tx)); ub = ((in_1)?(in_1):(&tz));
for (v = 0, b->n = ub->n = 0; v < n_vtx; ++v) {
@@ -11095,7 +11113,7 @@ uint32_t is_topo, uint32_t *max_drop_len)
// }
}
if(in_0) free(tx.a); if(in_1) free(tz.a);
if(in_0) {free(tx.a);} if(in_1) {free(tz.a);}
if (cnt > 0) usg_cleanup(g);
free(f);
// fprintf(stderr, "-[M::%s::] max_ext::%d, len_rat::%f\n", __func__, max_ext, len_rat);
@@ -11106,7 +11124,7 @@ uint32_t is_topo, uint32_t *max_drop_len, uint8_t *ff)
{
asg64_v tx = {0,0,0}, tz = {0,0,0}, *b = NULL, *ub = NULL;
uint32_t i, k, v, w, n_vtx = g->n<<1, nv, nw, kv, kw, /**trioF = (uint32_t)-1, ntrioF = (uint32_t)-1,**/ ol_max, ou_max, to_del, cnt = 0, mm_ol;
usg_arc_t *av, *aw, *ve, *we; uint64_t x, kocc[2]; uint8_t *f; CALLOC(f, g->n);
usg_arc_t *av = NULL, *aw = NULL, *ve = NULL, *we = NULL; uint64_t x, kocc[2]; uint8_t *f; CALLOC(f, g->n);
b = ((in_0)?(in_0):(&tx)); ub = ((in_1)?(in_1):(&tz));
for (v = 0, b->n = ub->n = 0; v < n_vtx; ++v) {
@@ -11225,7 +11243,7 @@ uint32_t is_topo, uint32_t *max_drop_len, uint8_t *ff)
// }
}
if(in_0) free(tx.a); if(in_1) free(tz.a);
if(in_0) {free(tx.a);} if(in_1) {free(tz.a);}
if (cnt > 0) usg_cleanup(g);
free(f);
// fprintf(stderr, "-[M::%s::] max_ext::%d, len_rat::%f\n", __func__, max_ext, len_rat);
@@ -11260,9 +11278,9 @@ uint64_t* nodeLen, uint64_t* baseLen, uint64_t *occ, asg64_v* b)
if(nodeLen) (*nodeLen) += g->a[v>>1].occ;
if(b) kv_push(uint64_t, *b, v);
if(occ) (*occ)++;
if(occ) {(*occ)++;}
///means reach the end of a unitig
if(kv!=1 && baseLen) (*baseLen) += g->a[v>>1].len;
if((kv!=1) && (baseLen)) {(*baseLen) += g->a[v>>1].len;}
if(kv==0) {
return_flag = END_TIPS; break;
}
@@ -11789,7 +11807,7 @@ void u2g_hybrid_extend(usg_t *ng, uint64_t* i_max_dist, asg64_v *in, asg64_v *ib
free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
if(n_pop) usg_cleanup(ng);
if(!in) free(tx.a); if(!ib) free(tb.a);
if(!in) {free(tx.a);} if(!ib) {free(tb.a);}
}
@@ -13748,7 +13766,7 @@ void u2g_hybrid_detan(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, asg64_v *
// fprintf(stderr, "-[M::%s::] ng->n::%u\n", __func__, (uint32_t)ng->n);
}
// prt_usg_t(uidx, ng, "ng_dbg");
if(!in) free(tx.a); if(!ib) free(tb.a);
if(!in) {free(tx.a);} if(!ib) {free(tb.a);}
}
void u2g_hybrid_aln(ul_resolve_t *uidx, usg_t *ng, asg64_v *ob, asg64_v *ub)
@@ -15384,7 +15402,7 @@ void u2g_hybrid_detan_iter(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint
// }
// prt_usg_t(uidx, ng, "ng1");
if(!in) free(tx.a); if(!ib) free(tb.a);
if(!in) {free(tx.a);} if(!ib) {free(tb.a);}
free(ng_occ); free(i_idx); free(ff); kv_destroy(b64); kv_destroy(ub64);
}
@@ -15826,7 +15844,7 @@ uint32_t is_topo, int32_t max_ext, float len_rat, uint8_t *f, asg64_v *in_0, asg
{
uint32_t i, k, v, w, nv, nw, kv, kw, ol_max, ou_max, to_del, cnt = 0, tip_v, mm_ol;
asg64_v tx = {0,0,0}, tz = {0,0,0}, *b = NULL, *ub = NULL; int32_t n_tip, r_tip;
uint64_t kocc[2]; usg_arc_t *av, *aw, *ve, *we;
uint64_t kocc[2]; usg_arc_t *av = NULL, *aw = NULL, *ve = NULL, *we = NULL;
b = ((in_0)?(in_0):(&tx)); ub = ((in_1)?(in_1):(&tz));
for (k = b->n = ub->n = 0; k < a_n; ++k) {
@@ -15893,10 +15911,10 @@ uint32_t is_topo, int32_t max_ext, float len_rat, uint8_t *f, asg64_v *in_0, asg
to_del = 1;
} else if (kw == 1) {
n_tip = usg_naive_topocut_aux(g, w^1, max_ext, f, b, ub);
if (n_tip < max_ext) to_del = 1; tip_v = w^1;
if (n_tip < max_ext) {to_del = 1;} tip_v = w^1;
} else if (kv == 1) {
n_tip = usg_naive_topocut_aux(g, v^1, max_ext, f, b, ub);
if (n_tip < max_ext) to_del = 1; tip_v = v^1;
if (n_tip < max_ext) {to_del = 1;} tip_v = v^1;
}
}
@@ -15910,7 +15928,7 @@ uint32_t is_topo, int32_t max_ext, float len_rat, uint8_t *f, asg64_v *in_0, asg
}
}
if(in_0) free(tx.a); if(in_1) free(tz.a);
if(in_0) {free(tx.a);} if(in_1) {free(tz.a);}
return cnt;
}
@@ -16719,7 +16737,7 @@ void renew_u2g_bg(ul_resolve_t *uidx)
for (v = 0, n_vtx = bg->bg->n_seq<<1; v < n_vtx; v++) {
if(bg->bg->seq[v>>1].del) continue;
av = asg_arc_a(bg->bg, v); nv = asg_arc_n(bg->bg, v);
if(!nv) continue; w = av[0].v;
if(!nv) {continue;} w = av[0].v;
v_occ[0] = v_occ[1] = w_occ[0] = w_occ[1] = (uint32_t)-1;
av = asg_arc_a(bg->bg, v); nv = asg_arc_n(bg->bg, v);
@@ -17845,7 +17863,7 @@ int usg_topocut_aux_del(ma_ug_t *ug, uint32_t v, int max_ext, asg64_v *b)
int32_t n_ext;
for (n_ext = 0; n_ext < max_ext && v != (uint32_t)-1; ) {
if (usg_topocut_aux_unambi1(ug->g, v^1) == (uint32_t)-1) break;
if(b) kv_push(uint64_t, *b, v); n_ext += ug->u.a[v>>1].n;
if(b) {kv_push(uint64_t, *b, v);} n_ext += ug->u.a[v>>1].n;
v = usg_topocut_aux_unambi1(ug->g, v);
}
return n_ext;
@@ -17894,7 +17912,7 @@ uint32_t min_node)
// fprintf(stderr, "[M::%s]\tStart\n", __func__);
ma_ug_t *ug = sl->ug; asg_t *g = sl->ug->g; uint32_t ol_max, ou_max, lnid;
uint32_t v, w, n_vtx = (g->n_seq<<1), nv, nw, i, k, z, kv, kw, bb, to_del, bn, tip;
asg64_v tx = {0,0,0}, *b = NULL; asg_arc_t *av, *aw, *ve, *we; ma_utg_t *u;
asg64_v tx = {0,0,0}, *b = NULL; asg_arc_t *av = NULL, *aw = NULL, *ve = NULL, *we = NULL; ma_utg_t *u;
uint32_t trioF = (uint32_t)-1, ntrioF = (uint32_t)-1, mm_ol, mm_ou, cnt = 0, del_v, del_w;
if(in) b = in;
hic.cpp
+1 -1
View File
@@ -6569,7 +6569,7 @@ min_cut_t* m, hc_links* link, G_partition* x)
if(res->full_bub == 0)
{
res->a.n = 0;
uint32_t v, u = 0, uv, k_n, pre_n = x->n;
uint32_t v, u = 0, uv = UINT32_MAX, k_n, pre_n = x->n;
hc_linkeage* t = NULL;
x->n--;
for (i = 0; i < n; i++)
horder.cpp
+2 -2
View File
@@ -32,7 +32,7 @@ KRADIX_SORT_INIT(osg, osg_arc_t, osg_arc_key, member_size(osg_arc_t, u))
#define BREAK_CUTOFF 0.1
#define BREAK_BOUNDARY 0.015
void reduce_hamming_error_adv(ma_ug_t *iug, asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub);
int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub, uint32_t max_ext);
typedef struct {
uint64_t ruid;
@@ -3936,7 +3936,7 @@ asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, kv_u_trans_t *ref, ug_opt_t *opt,
// output_hic_rtg(i_ug, h->r_g, opt, asm_opt.output_file_name);
// reduce_hamming_error(h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz);
reduce_hamming_error_adv(NULL, h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz, opt->ruIndex, NULL);
reduce_hamming_error_adv(NULL, h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz, opt->ruIndex, NULL, (asm_opt.max_short_tip*2));
/**
scaffold_hap(h, t_idx, opt, round, asm_opt.output_file_name, FATHER);
scaffold_hap(h, t_idx, opt, round, asm_opt.output_file_name, MOTHER);
htab.cpp
+120 -18
View File
@@ -385,7 +385,7 @@ ha_pt_t *ha_pt_gen(ha_ct_t *ct, int n_thread, int is_l)
ha_ct_destroy_bf(ct);
CALLOC(pt, 1);
pt->k = ct->k, pt->pre = ct->pre, pt->tot = ct->tot;
CALLOC(pt->h, 1<<pt->pre);
CALLOC(pt->h, (((uint64_t)1)<<pt->pre));
for (i = 0; i < 1<<pt->pre; ++i) {
pt->h[i].h = yak_pt_init();
yak_pt_resize(pt->h[i].h, kh_size(ct->h[i].h));
@@ -422,7 +422,7 @@ ha_pt_t *ha_pt_gen_count(ha_ct_t *ct, int n_thread)
ha_ct_destroy_bf(ct);
CALLOC(pt, 1);
pt->k = ct->k, pt->pre = ct->pre, pt->tot = ct->tot;
CALLOC(pt->h, 1<<pt->pre);
CALLOC(pt->h, (((uint64_t)1)<<pt->pre));
for (i = 0; i < 1<<pt->pre; ++i) {
pt->h[i].h = yak_pt_init();
yak_pt_resize(pt->h[i].h, kh_size(ct->h[i].h));
@@ -578,7 +578,7 @@ KSEQ_INIT(gzFile, gzread)
typedef struct { // global data structure for kt_pipeline()
const yak_copt_t *opt;
const void *flt_tab;
int flag, create_new, is_store, uq;
int flag, create_new, is_store, uq, ifq;
uint64_t n_mz, n_seq; ///number of total reads
kseq_t *ks;
UC_Read ucr;
@@ -761,7 +761,7 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for
while ((ret = kseq_read(p->ks)) >= 0) {\
int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen);\
if((l <= 0) || (l < asm_opt.rl_cut)) continue;\
if((asm_opt.is_sc) && (asm_opt.sc_cut > 0) && (!flt_quals(p->ks->qual.s+p->opt->adaLen, l, 33, asm_opt.sc_cut))) continue;\
if((p->ifq) && (asm_opt.sc_cut > 0) && (!flt_quals(p->ks->qual.s+p->opt->adaLen, l, 33, asm_opt.sc_cut))) continue;\
if (p->n_seq >= 1<<28) {\
fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);\
exit(1);\
@@ -779,7 +779,7 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for
++n_N;\
ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N);\
memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);\
if(p->rs_out->rsc) {\
if(p->ifq) {\
ha_compress_qual(Get_QUAL(*p->rs_out, p->n_seq), p->ks->qual.s+p->opt->adaLen, l, sc_bn, 33);\
/**print_fastq(NULL, p->ks->name.s, p->ks->seq.s, p->ks->qual.s, (1<<sc_bn), 33);**/\
/**if(l <= 1000000) {\
@@ -836,7 +836,7 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for
uint32_t j;\
/**s->n_seq is how many reads at this buffer**/\
/**s->mz && s->mz_buf are lists of minimzer vectors**/\
CALLOC(s->mz, s->n_seq), CALLOC(s->mz_buf, p->opt->n_thread), CALLOC(s->mt, p->opt->n_thread);\
CALLOC(s->mz, s->n_seq); CALLOC(s->mz_buf, p->opt->n_thread); CALLOC(s->mt, p->opt->n_thread);\
/**calculate minimzers for each read, each read corresponds to one thread**/\
kt_for(p->opt->n_thread, sf##_worker_for_mz, s, s->n_seq);\
for (i = 0; i < p->opt->n_thread; ++i) free(s->mt[i].a), free(s->mz_buf[i].a);\
@@ -927,7 +927,7 @@ void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs)
exit(1);
}
static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt_t *p0, ha_ct_t *c0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int64_t *n_seq)
static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt_t *p0, ha_ct_t *c0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int64_t *n_seq, uint8_t ifq)
{
///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN|HAF_CREATE_NEW
int read_rs = (rs && (flag & HAF_RS_READ));
@@ -935,7 +935,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
pl_data_t pl;
gzFile fp = 0;
memset(&pl, 0, sizeof(pl_data_t));
pl.n_seq = *n_seq;
pl.n_seq = *n_seq; pl.ifq = ifq;
if(ug_rs) {
pl.us_in = us;
} else if (read_rs) {
@@ -1019,9 +1019,26 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_o, int flag, int HPC, int k, int w, h
}**/
///asm_opt->num_reads is the number of fastq files
for (i = 0; i < (us?1:asm_o->num_reads); ++i){
h = yak_count(&opt, asm_o->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq);
h = yak_count(&opt, asm_o->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq, asm_opt.is_sc);
if(h) n_bs += h->bs;
}
if((rs) && (flag & HAF_RS_WRITE_LEN) && (asm_opt.is_sc)) {
rs->tqn = rs->total_reads; rs->tr[0] = rs->total_reads_bases;
}
if((asm_o->hf) && (!us)) {
for (i = 0; i < asm_o->hf->n; ++i) {
h = yak_count(&opt, asm_o->hf->a[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq, 0);
if(h) n_bs += h->bs;
}
}
if((rs) && (flag & HAF_RS_WRITE_LEN) && (asm_opt.is_sc)) {
rs->tr[1] = rs->total_reads_bases - rs->tr[0];
}
// fprintf(stderr, "[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]);
if(h) h->bs = n_bs;
if (h && opt.bf_shift > 0)
ha_ct_destroy_bf(h);
@@ -1141,7 +1158,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i
if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ;
ha_ct_t *h;
h = ha_count(asm_opt, HAF_COUNT_ALL|ex_flag|((read_from_store)?(HAF_RS_READ):(HAF_RS_WRITE_LEN)), !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, NULL, rs, NULL, 1, NULL, 0);
if((asm_opt->flag & HA_F_VERBOSE_GFA))
if((asm_opt->flag & HA_F_VERBOSE_GFA) || (asm_opt->restart))
{
write_ct_index((void*)h, asm_opt->output_file_name);
// load_ct_index(&ha_ct_table, asm_opt->output_file_name);
@@ -1154,10 +1171,20 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i
{
ha_ct_hist(h, cnt, asm_opt->thread_num);
peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, asm_opt->hg_size>0?(h->bs/asm_opt->hg_size):(-1), cnt, &peak_het);
///r850
fprintf(stderr, "[M::%s::auto] inferred peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
if((asm_opt->het_cov_ss > 0) || (asm_opt->hmo_cov_ss > 0)) {
fprintf(stderr, "[M::%s::user] override requested peak_hom: %ld; peak_het: %ld\n", __func__, asm_opt->hmo_cov_ss, asm_opt->het_cov_ss);
if(asm_opt->hmo_cov_ss > 0) peak_hom = asm_opt->hmo_cov_ss;
if(asm_opt->het_cov_ss > 0) peak_het = asm_opt->het_cov_ss;
}
if (hom_cov) *hom_cov = peak_hom;
if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
if (peak_hom > 0) fprintf(stderr, "[M::%s::final] using peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
///in default, asm_opt->high_factor = 5.0
///r833
cutoff = (int)(peak_hom * asm_opt->high_factor);
if(cutoff < asm_opt->hf_cutoff) cutoff = asm_opt->hf_cutoff;
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
}
ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num);
@@ -1252,9 +1279,16 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
ha_ct_hist(ct, cnt, asm_opt->thread_num);
fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]);
peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, asm_opt->hg_size>0?(ct->bs/asm_opt->hg_size):(-1), cnt, &peak_het);
///r850
fprintf(stderr, "[M::%s::auto] inferred peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
if((asm_opt->het_cov_ss > 0) || (asm_opt->hmo_cov_ss > 0)) {
fprintf(stderr, "[M::%s::user] override requested peak_hom: %ld; peak_het: %ld\n", __func__, asm_opt->hmo_cov_ss, asm_opt->het_cov_ss);
if(asm_opt->hmo_cov_ss > 0) peak_hom = asm_opt->hmo_cov_ss;
if(asm_opt->het_cov_ss > 0) peak_het = asm_opt->het_cov_ss;
}
if (hom_cov) *hom_cov = peak_hom;
if (het_cov) *het_cov = peak_het;
if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
if (peak_hom > 0) fprintf(stderr, "[M::%s::final] using peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
///here ha_ct_shrink is mostly used to remove k-mer appearing only 1 time
if (flt_tab == 0) {
int cutoff = (int)(peak_hom * asm_opt->high_factor);
@@ -1367,7 +1401,8 @@ int load_ct_index(void **i_ct_idx, 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)
{
char* gfa_name = (char*)malloc(strlen(file_name)+64);
sprintf(gfa_name, "%s.pt_flt", file_name);
if(r) sprintf(gfa_name, "%s.pt_flt", file_name);
else sprintf(gfa_name, "%s.pt_flt.bin", file_name);
FILE* fp = fopen(gfa_name, "w");
if (!fp) {
free(gfa_name);
@@ -1406,7 +1441,7 @@ int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t*
fwrite(&opt->het_cov, sizeof(opt->het_cov), 1, fp);
fwrite(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp);
if(r) {
write_All_reads(r, gfa_name);
sprintf(gfa_name, "%s.pt_flt.paf.bin", file_name);
@@ -1422,6 +1457,7 @@ 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);
}
}
fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
free(gfa_name);
@@ -1432,7 +1468,9 @@ int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t*
int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_opt_t* opt, char* file_name)
{
char* gfa_name = (char*)malloc(strlen(file_name)+64);
sprintf(gfa_name, "%s.pt_flt", file_name);
if(r) sprintf(gfa_name, "%s.pt_flt", file_name);
else sprintf(gfa_name, "%s.pt_flt.bin", file_name);
// fprintf(stderr, "[M::%s]\tgfa_name::%s\n", __func__, gfa_name);
FILE* fp = fopen(gfa_name, "r");
if (!fp) {
free(gfa_name);
@@ -1512,9 +1550,8 @@ int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_op
// fclose(fp);
if(!load_All_reads(r, gfa_name))
{
if(r) {
if(!load_All_reads(r, gfa_name)) {
free(gfa_name);
return 0;
}
@@ -1545,6 +1582,7 @@ 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);
}
}
fclose(fp);
fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__);
@@ -1553,6 +1591,70 @@ int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_op
return 1;
}
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)
{
char* gfa_name = (char*)malloc(strlen(file_name)+64);
sprintf(gfa_name, "%s.ad", file_name);
if(is_w) {
write_pt_index(*flt_tab, *ha_idx, NULL, opt, gfa_name);
} else {
load_pt_index(flt_tab, ha_idx, NULL, opt, gfa_name);
}
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)
{
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");
if (!fp) {
free(gfa_name);
return 0;
}
f_flag += fread(&rr0, sizeof(rr0), 1, fp);
f_flag += fread(&tot_rr0, sizeof(tot_rr0), 1, fp);
fclose(fp);
if(rr0 != rr || tot_rr0 != tot_rr) {
free(gfa_name);
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)) {
free(gfa_name);
return 0;
}
} else {
sprintf(gfa_name, "%s.r%lu", file_name, rr);
write_pt_index(*r_flt_tab, *r_ha_idx, r, opt, gfa_name);
sprintf(gfa_name, "%s.r%lu.ht.bin", file_name, rr);
fp = fopen(gfa_name, "w");
if (!fp) {
free(gfa_name);
return 0;
}
fwrite(&rr, sizeof(rr), 1, fp);
fwrite(&tot_rr, sizeof(tot_rr), 1, fp);
fclose(fp);
}
free(gfa_name);
return 1;
}
int uidx_write(void *flt_tab, ha_pt_t *ha_idx, char* file_name, ma_ug_t *ug)
htab.h
+4 -1
View File
@@ -88,11 +88,13 @@ 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);
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);
ha_abuf_t *ha_abuf_init_buf(void *km);
ha_abufl_t *ha_abufl_init_buf(void *km);
@@ -117,6 +119,7 @@ double yak_cpu_usage(void);
void ha_triobin(const hifiasm_opt_t *opt);
uint32_t test_yak_binning(char* fn, char *cmd);
uint32_t *ha_polybin_list(const hifiasm_opt_t *opt);
uint32_t *ha_charbin_list(const hifiasm_opt_t *opt, uint8_t **idx, uint32_t *idx_n);
void mz1_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique, void *km);
void mz2_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mzl_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique, void *km);
@@ -166,7 +169,7 @@ static inline uint64_t yak_hash_long(uint64_t x[4])
return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]);
}
#define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr))))
#define CALLOC(ptr, len) ((ptr) = ((((len)*sizeof(*(ptr))) <= 9223372036854775807)?((__typeof__(ptr))calloc((len), sizeof(*(ptr)))):(NULL)))
#define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr))))
#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr))))
#define MEMCPY(dest, src, len) (memcpy((dest), (src), (len) * sizeof(*(src))))
inter.cpp
+22 -22
View File
@@ -1381,8 +1381,8 @@ const int64_t qlen, const int64_t rlen, int32_t *r_qs, int32_t *r_qe, int32_t *r
re = rlen - 1; qe += rtail;
}
if(r_qs) (*r_qs) = qs; if(r_qe) (*r_qe) = qe + 1;
if(r_rs) (*r_rs) = rs; if(r_re) (*r_re) = re + 1;
if(r_qs) {(*r_qs) = qs;} if(r_qe) {(*r_qe) = qe + 1;}
if(r_rs) {(*r_rs) = rs;} if(r_re) {(*r_re) = re + 1;}
if(rev) {
if(r_rs) (*r_rs) = rlen - re - 1;
if(r_re) (*r_re) = rlen - rs;
@@ -1413,8 +1413,8 @@ int32_t *r_qs, int32_t *r_qe, int32_t *r_rs, int32_t *r_re)
re = rlen - 1; qe += rtail;
}
if(r_qs) (*r_qs) = qs; if(r_qe) (*r_qe) = qe + 1;
if(r_rs) (*r_rs) = rs; if(r_re) (*r_re) = re + 1;
if(r_qs) {(*r_qs) = qs;} if(r_qe) {(*r_qe) = qe + 1;}
if(r_rs) {(*r_rs) = rs;} if(r_re) {(*r_re) = re + 1;}
if(ri->v&1) {
if(r_rs) (*r_rs) = rlen - re - 1;
if(r_re) (*r_re) = rlen - rs;
@@ -4500,7 +4500,7 @@ int64_t g_adjacent_dis(const asg_t *g, uint32_t v, uint32_t w)
void get_r_offset(ma_ug_t *ug, mg_lchain_t *x, int64_t *rs, int64_t *re, int64_t *qs, int64_t *qe)
{
if(qs) *qs = x->qs; if(qe) *qe = x->qe;
if(qs) {*qs = x->qs;} if(qe) {*qe = x->qe;}
if(!(x->score&1)) {
if(rs) *rs = x->rs + x->off;
if(re) *re = x->re + x->off;
@@ -4512,7 +4512,7 @@ void get_r_offset(ma_ug_t *ug, mg_lchain_t *x, int64_t *rs, int64_t *re, int64_t
void get_u_offset(ma_ug_t *ug, mg_lchain_t *x, int64_t *rs, int64_t *re, int64_t *qs, int64_t *qe)
{
if(qs) *qs = x->qs; if(qe) *qe = x->qe;
if(qs) {*qs = x->qs;} if(qe) {*qe = x->qe;}
if(!(x->v&1)) {
if(rs) *rs = x->rs + x->off;
if(re) *re = x->re + x->off;
@@ -6730,13 +6730,13 @@ int64_t comput_sc_partial_cigar(int64_t sc, int64_t ol, double err_sc_r, overlap
// pk = (*wi); pe = (*werr);
if(k < 0) {
k = 0; err = 0;
if(wi) (*wi) = k; if(werr) (*werr) = err;
if(wi) {(*wi) = k;} if(werr) {(*werr) = err;}
// assert(e <= z->w_list.a[0].x_start);
} else {
if(z->w_list.a[k].y_end != -1) {
err -= z->w_list.a[k].error;
}
if(wi) (*wi) = k; if(werr) (*werr) = err;
if(wi) {(*wi) = k;} if(werr) {(*werr) = err;}
// if(!(err >= 0 && k >= 0 && k < wn && z->w_list.a[k].x_start < e && z->w_list.a[k].x_end + 1 >= e)){
// fprintf(stderr, "[M::%s] ol::%ld, e::%ld, z::[%u, %u], k::%ld, wn::%ld, w::[%d, %d], err::%ld\n", __func__,
@@ -8817,7 +8817,7 @@ double e_rate, kv_rtrace_t *trace, uint32_t rechain_w, ul_ov_t *res, ul_ov_t *rr
for (i = k = 0; i < nv; i++) {
if(res[i].qn == (uint32_t)-1) continue;
res[k] = res[i];
if(res[k].sec&mm) res[k].sec -= mm; k++;
if(res[k].sec&mm) {res[k].sec -= mm;} k++;
// res[k] = res[i];
// id[k] = res[k].sec;
// if(id[k]&mm) id[k] -= mm;
@@ -10099,7 +10099,7 @@ uint32_t test_het_aln(ma_ug_t *ug, uint64_t rid, u_trans_t *a, uint64_t a_n, st_
uint32_t is_mmhom_node(uint64_t *ca, ma_utg_t *u, asg_t *sg, uint64_t cov_bd, double cut_rate)
{
if(cut_rate < 0) cut_rate = 0; if(cut_rate > 1.0) cut_rate = 1.0;
if(cut_rate < 0) {cut_rate = 0;} if(cut_rate > 1.0) {cut_rate = 1.0;}
uint64_t k, a, na, a_cut = u->n*cut_rate, na_cut = u->n*(1.0-cut_rate);
for (k = a = na = 0; k < u->n; k++) {
if(ca[k] > (cov_bd*((uint64_t)sg->seq[u->a[k]>>33].len))) {
@@ -13637,8 +13637,8 @@ void gen_end_coord(ul_ov_t *z, int64_t qlen, int64_t tlen, int64_t *r_qs, int64_
te = tlen; qe += ttail;
}
if(r_qs) (*r_qs) = qs; if(r_qe) (*r_qe) = qe;
if(r_ts) (*r_ts) = ts; if(r_te) (*r_te) = te;
if(r_qs) {(*r_qs) = qs;} if(r_qe) {(*r_qe) = qe;}
if(r_ts) {(*r_ts) = ts;} if(r_te) {(*r_te) = te;}
if(z->rev) {
if(r_ts) (*r_ts) = tlen - te;
if(r_te) (*r_te) = tlen - ts;
@@ -13822,8 +13822,8 @@ void extend_end_coord(mg_lchain_t *li, ul_ov_t *ui, const int64_t qlen, const in
re = rlen; qe += rtail;
}
if(r_qs) (*r_qs) = qs; if(r_qe) (*r_qe) = qe;
if(r_rs) (*r_rs) = rs; if(r_re) (*r_re) = re;
if(r_qs) {(*r_qs) = qs;} if(r_qe) {(*r_qe) = qe;}
if(r_rs) {(*r_rs) = rs;} if(r_re) {(*r_re) = re;}
if(rev) {
if(r_rs) (*r_rs) = rlen - re;
if(r_re) (*r_re) = rlen - rs;
@@ -14327,7 +14327,7 @@ uint64_t detect_mul_way, float len_dif)
int64_t hc_chain_backtrack(int64_t n, const int64_t *f, const uint64_t *p, uint64_t *srt, uint64_t *u, uint64_t *v,
int64_t *n_u_, int64_t *n_v_)
{
if(n_u_) *n_u_ = 0; if(n_v_) *n_v_ = 0;
if(n_u_) {*n_u_ = 0;} if(n_v_) {*n_v_ = 0;}
int64_t i, k, n_v, n_srt, n_v0, n_u, sc;
if (n == 0) return 0;
// v[] keeps the peak score up to i; f[] is the score ending at i, not always the peak
@@ -14351,7 +14351,7 @@ int64_t *n_u_, int64_t *n_v_)
u[n_u++] = (((uint64_t)sc)<<32) | ((uint64_t)(n_v-n_v0));
}
if(n_u_) *n_u_ = n_u; if(n_v_) *n_v_ = n_v;
if(n_u_) {*n_u_ = n_u;} if(n_v_) {*n_v_ = n_v;}
return n_u;
}
@@ -15532,7 +15532,7 @@ void update_rovlp_chain_qse(ma_ug_t *ug, int64_t sidx, int64_t eidx, mg_lchain_t
for (i = sidx+1; i < eidx; i++) {
get_r_offset(ug, &(a[i]), &rs, &re, NULL, NULL);
a[i].qe = cal_qext_coor(left_r[1], re, left_q[1], left_q[1] + re - left_r[1], re);
if(a[i].qe < 0) a[i].qe = 0; if(a[i].qe > qlen) a[i].qe = qlen;
if(a[i].qe < 0) {a[i].qe = 0;} if(a[i].qe > qlen) {a[i].qe = qlen;}
a[i].qs = cal_qext_coor(left_r[0], (a[i].qe<=left_q[1])?re:left_r[1],
left_q[0], (a[i].qe<=left_q[1])?a[i].qe:left_q[1], rs);
assert(a[i].qs >= 0 && a[i].qs <= qlen);
@@ -15553,7 +15553,7 @@ void update_rovlp_chain_qse(ma_ug_t *ug, int64_t sidx, int64_t eidx, mg_lchain_t
a[i].qe = cal_qext_coor(right_r[0], right_r[1], right_q[0], right_q[1], re);
assert(a[i].qe >= 0 && a[i].qe <= qlen);
a[i].qs = cal_qext_coor(rs, right_r[0], right_q[0]-(right_r[0]-rs), right_q[0], rs);
if(a[i].qs < 0) a[i].qs = 0; if(a[i].qs > qlen) a[i].qs = qlen;
if(a[i].qs < 0) {a[i].qs = 0;} if(a[i].qs > qlen) {a[i].qs = qlen;}
if(a[i].qs > a[i].qe) {
tt = a[i].qs; a[i].qs = a[i].qe; a[i].qe = tt;
}
@@ -16345,8 +16345,8 @@ void ctg_rg2ug_gen(ma_utg_t *r_cl, kv_ul_ov_t *u_cl, const ul_idx_t *uref, uint6
for (i = l; i < k; i++) {///could be merged
m = get_unique_rctg_aln(uref, r_cl->a[i]>>32, i, ls, r_cl->len, &kp, NULL, ((uint32_t)-1));
assert(m == 1); assert(kp.tn == lp.tn); assert(kp.qn == lp.qn + i - l); assert(kp.sec == lp.sec + i - l);
if(kp.qs < lp.qs) lp.qs = kp.qs; if(kp.qe > lp.qe) lp.qe = kp.qe;
if(kp.ts < lp.ts) lp.ts = kp.ts; if(kp.te > lp.te) lp.te = kp.te;
if(kp.qs < lp.qs) {lp.qs = kp.qs;} if(kp.qe > lp.qe) {lp.qe = kp.qe;}
if(kp.ts < lp.ts) {lp.ts = kp.ts;} if(kp.te > lp.te) {lp.te = kp.te;}
ls += (uint32_t)(r_cl->a[i]);
// if(id == 123) {
@@ -18590,7 +18590,7 @@ void collect_pp_ovlps(ul_ov_t *a, uint64_t a_n, kv_ul_ov_t *res, const ug_opt_t
if(((*mqe) == ((uint64_t)-1)) || ((*mqe) < a[i].qe)) (*mqe) = a[i].qe;
}
if(st) (*mqs) = st->qs; if(et) (*mqe) = et->qe;
if(st) {(*mqs) = st->qs;} if(et) {(*mqe) = et->qe;}
}
uint64_t gen_cns_chain_linear_hard(ul_ov_t *a, int64_t a_n, const asg_t *rg, int64_t qlen, int64_t bw, double diff_thre, Chain_Data* dp,
@@ -21768,7 +21768,7 @@ void asg_arc_push_contain_trans(asg_t *g, ma_hit_t_alloc* ov, int64_t min_ovlp,
nv = asg_arc_n(g, v); av = asg_arc_a(g, v); avi = g->idx[v]>>32;
while(nv) {
for (i = nc = cc = 0; i < nv; ++i) {
if(av[i].del) continue; w = av[i].v;
if(av[i].del) {continue;} w = av[i].v;
if(!(is_contain_r((*ri), (w>>1)))) continue;///new arcs must be bridged by contained reads
assert(!(g->seq[w>>1].del));
nw = asg_arc_n(g, w); aw = asg_arc_a(g, w); awi = g->idx[w]>>32;
rcut.cpp
+15 -6
View File
@@ -806,8 +806,19 @@ uint32_t mc_edges_symm(mc_match_t *ma)
static void normalize_mb_edge(mb_edge_t *a, mb_edge_t *b)
{
if(a->w >= b->w)
{
uint8_t f = 0;
if(a->w[0] > b->w[0]) {
f = 1;
} else if(a->w[0] == b->w[0] && a->w[1] > b->w[1]) {
f = 1;
} else if(a->w[0] == b->w[0] && a->w[1] == b->w[1] && a->w[2] > b->w[2]) {
f = 1;
} else if(a->w[0] == b->w[0] && a->w[1] == b->w[1] && a->w[2] == b->w[2] && a->w[3] > b->w[3]) {
f = 1;
}
// if(a->w >= b->w) {
if(f) {
b->x = (uint32_t)a->x;
b->x <<= 32;
b->x |= (a->x>>32);
@@ -815,9 +826,7 @@ static void normalize_mb_edge(mb_edge_t *a, mb_edge_t *b)
b->w[3] = a->w[3];
b->w[1] = a->w[2];
b->w[2] = a->w[1];
}
else
{
} else {
a->x = (uint32_t)b->x;
a->x <<= 32;
a->x |= (b->x>>32);
@@ -3346,7 +3355,7 @@ mc_clus_t *init_mc_clus_t(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub, uin
mc_clus_t *p; CALLOC(p, 1);
p->bub = bub; p->opt = opt; p->mg = mg; p->n = bub->ug->g->n_seq;
CALLOC(p->lock, p->n);
if(n_thread > 64) n_thread = 64; p->n_thread = n_thread;
if(n_thread > 64) {n_thread = 64;} p->n_thread = n_thread;
CALLOC(p->aux, p->n_thread);
uint32_t k, ss = (p->n>>3)+(!!(p->n&7));
for (k = 0; k < p->n_thread; k++) {