mirror of
https://github.com/chhylp123/hifiasm.git
synced 2026-09-28 13:38:11 +08:00
for debugging
This commit is contained in:
chhylp123
+21
-9
@@ -12,7 +12,7 @@
|
||||
#include "kthread.h"
|
||||
|
||||
void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres,
|
||||
int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar);
|
||||
int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct);
|
||||
void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list);
|
||||
|
||||
All_reads R_INF;
|
||||
@@ -561,7 +561,7 @@ static void worker_ovec(void *data, long i, int tid)
|
||||
int fully_cov, abnormal;
|
||||
|
||||
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist,
|
||||
0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region));
|
||||
0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
|
||||
|
||||
clear_Cigar_record(&b->cigar1);
|
||||
clear_Round2_alignment(&b->round2);
|
||||
@@ -606,7 +606,6 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
|
||||
uint64_t k, queryNameLen;
|
||||
for (k = 0; k < R_INF_FLAG.query_num; k++)
|
||||
{
|
||||
|
||||
queryNameLen = strlen(R_INF_FLAG.read_name[k]);
|
||||
if (queryNameLen != Get_NAME_LENGTH((R_INF),i)) continue;
|
||||
if (memcmp(R_INF_FLAG.read_name[k], Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
|
||||
@@ -617,10 +616,10 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
|
||||
|
||||
if(k < R_INF_FLAG.query_num)
|
||||
{
|
||||
int fully_cov, abnormal;
|
||||
int fully_cov, abnormal, q_idx = k;
|
||||
|
||||
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist,
|
||||
0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region));
|
||||
0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), &(R_INF_FLAG.candidate_count[q_idx]));
|
||||
|
||||
clear_Cigar_record(&b->cigar1);
|
||||
clear_Round2_alignment(&b->round2);
|
||||
@@ -649,6 +648,18 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
|
||||
fprintf(R_INF_FLAG.fp, "\n>%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i));
|
||||
fprintf(R_INF_FLAG.fp, "%d-th round, len: %lu, hom_cov: %d, max_n_chain: %d\n",
|
||||
asm_opt.number_of_round, Get_READ_LENGTH(R_INF, i), asm_opt.hom_cov, asm_opt.max_n_chain);
|
||||
|
||||
fprintf(R_INF_FLAG.fp, "***************************k-mer counts (%d)***************************\n", (int)(R_INF_FLAG.candidate_count[q_idx].a.n));
|
||||
|
||||
sort_kvec_t_u64_warp(&(R_INF_FLAG.candidate_count[q_idx]), 0);
|
||||
for (k = 0; k < R_INF_FLAG.candidate_count[q_idx].a.n; k++)
|
||||
{
|
||||
fprintf(R_INF_FLAG.fp, "[%lu] Count(%u): %lu, filtered: %lu\n", k,
|
||||
(uint32_t)R_INF_FLAG.candidate_count[q_idx].a.a[k], R_INF_FLAG.candidate_count[q_idx].a.a[k]>>33,
|
||||
(R_INF_FLAG.candidate_count[q_idx].a.a[k]>>32)&(uint64_t)1);
|
||||
}
|
||||
|
||||
|
||||
fprintf(R_INF_FLAG.fp, "***************************forward ovlp***************************\n");
|
||||
for (k = 0; k < b->olist.length; k++)
|
||||
{
|
||||
@@ -692,6 +703,7 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
|
||||
// }
|
||||
|
||||
fprintf(R_INF_FLAG.fp, "R_INF.trio_flag[%ld]: %u\n", i, R_INF.trio_flag[i]);
|
||||
|
||||
|
||||
pthread_mutex_unlock(&R_INF_FLAG.OutputMutex);
|
||||
}
|
||||
@@ -867,7 +879,7 @@ void ha_overlap_and_correct(int round)
|
||||
else
|
||||
kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads);
|
||||
|
||||
if (r_out) write_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);
|
||||
ha_pt_destroy(ha_idx);
|
||||
ha_idx = NULL;
|
||||
|
||||
@@ -1298,7 +1310,7 @@ static void worker_ov_final(void *data, long i, int tid)
|
||||
|
||||
//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
|
||||
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 0.001,
|
||||
asm_opt.max_n_chain, 0, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region));
|
||||
asm_opt.max_n_chain, 0, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
|
||||
|
||||
overlap_region_sort_y_id(b->olist.list, b->olist.length);
|
||||
ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length);
|
||||
@@ -1364,7 +1376,7 @@ static void worker_ov_final_high_het(void *data, long i, int tid)
|
||||
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
|
||||
|
||||
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, HIGH_HET_ERROR_RATE,
|
||||
asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region));
|
||||
asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
|
||||
|
||||
overlap_region_sort_y_id(b->olist.list, b->olist.length);
|
||||
ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length);
|
||||
@@ -1631,7 +1643,7 @@ int ha_assemble(void)
|
||||
}
|
||||
if (!ovlp_loaded) {
|
||||
ha_flt_tab = ha_idx = NULL;
|
||||
if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
|
||||
if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name), load_ct_index(&ha_ct_table, asm_opt.output_file_name);
|
||||
|
||||
// construct hash table for high occurrence k-mers
|
||||
if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL)
|
||||
|
||||
+1
-1
@@ -3,7 +3,7 @@
|
||||
|
||||
#include <pthread.h>
|
||||
|
||||
#define HA_VERSION "0.12-r304"
|
||||
#define HA_VERSION "0.12-r305"
|
||||
|
||||
#define VERBOSE 0
|
||||
|
||||
|
||||
+25
-9
@@ -2107,20 +2107,30 @@ static inline int asg_is_single_edge(const asg_t *g, uint32_t v, uint32_t start_
|
||||
return nv;
|
||||
}
|
||||
|
||||
void debug_info_of_specfic_node(const char* name, asg_t *g, char* command)
|
||||
void debug_info_of_specfic_node(char* name, asg_t *g, R_to_U* ruIndex, char* command)
|
||||
{
|
||||
fprintf(stderr, "\n\n\n");
|
||||
uint32_t v, n_vtx = g->n_seq * 2, queryLen = strlen(name), flag = 0;
|
||||
uint32_t v, n_vtx = g->n_seq * 2, queryLen = strlen(name), flag = 0, contain_rId, is_Unitig;
|
||||
for (v = 0; v < n_vtx; ++v)
|
||||
{
|
||||
if(queryLen == Get_NAME_LENGTH(R_INF, (v>>1)) && memcmp(name, Get_NAME(R_INF, (v>>1)), Get_NAME_LENGTH(R_INF, (v>>1))) == 0)
|
||||
{
|
||||
if(flag == 0) fprintf(stderr, "\nafter %s\n", command);
|
||||
fprintf(stderr, "****************graph ref_read: %.*s, dir: %u****************\n",
|
||||
(int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), v&1);
|
||||
fprintf(stderr, "****************graph ref_read: %.*s, id: %u, dir: %u****************\n",
|
||||
(int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), v>>1, v&1);
|
||||
if(g->seq[v>>1].del)
|
||||
{
|
||||
fprintf(stderr, "read has already been deleted.\n");
|
||||
get_R_to_U(ruIndex, (v>>1), &contain_rId, &is_Unitig);
|
||||
if(contain_rId != (uint32_t)-1 && is_Unitig != 1)
|
||||
{
|
||||
fprintf(stderr, "read is deleted as a contained read by: %.*s\n contain_rId: %u, del: %u\n",
|
||||
(int)Get_NAME_LENGTH(R_INF, contain_rId), Get_NAME(R_INF, contain_rId), contain_rId, g->seq[contain_rId].del);
|
||||
}
|
||||
else
|
||||
{
|
||||
fprintf(stderr, "read has already been deleted.\n");
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
@@ -9757,12 +9767,13 @@ ma_hit_t_alloc* reverse_sources, int id, const char* command)
|
||||
for (j = 0; j < reverse_sources[i].length; j++)
|
||||
{
|
||||
tn = Get_tn(reverse_sources[i].buffer[j]);
|
||||
fprintf(stderr, "target: %.*s, qs: %u, qe: %u, ts: %u, te: %u, del: %u\n",
|
||||
fprintf(stderr, "target: %.*s, qs: %u, qe: %u, ts: %u, te: %u, rev: %u, del: %u\n",
|
||||
(int)Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn),
|
||||
Get_qs(reverse_sources[i].buffer[j]),
|
||||
Get_qe(reverse_sources[i].buffer[j]),
|
||||
Get_ts(reverse_sources[i].buffer[j]),
|
||||
Get_te(reverse_sources[i].buffer[j]),
|
||||
reverse_sources[i].buffer[j].rev,
|
||||
(uint32_t)sources[i].buffer[j].del);
|
||||
}
|
||||
|
||||
@@ -26919,6 +26930,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
|
||||
///it seems we do not need ma_hit_cut & ma_hit_flt
|
||||
ma_hit_sub(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut);
|
||||
detect_chimeric_reads(sources, n_read, readLen, coverage_cut, asm_opt.max_ov_diff_final * 2.0);
|
||||
|
||||
ma_hit_cut(sources, n_read, readLen, mini_overlap_length, &coverage_cut);
|
||||
///print_binned_reads(sources, n_read, coverage_cut);
|
||||
ma_hit_flt(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
|
||||
@@ -26926,8 +26938,10 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
|
||||
///fix_binned_reads(sources, n_read, coverage_cut);
|
||||
///just need to deal with trio here
|
||||
ma_hit_contained_advance(sources, n_read, coverage_cut, ruIndex, max_hang_length, mini_overlap_length);
|
||||
|
||||
sg = ma_sg_gen(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
|
||||
|
||||
///debug_info_of_specfic_node((char*)"m64062_190804_172951/130483063/ccs", sg, ruIndex, (char*)"sbsbsb");
|
||||
|
||||
asg_arc_del_trans(sg, gap_fuzz);
|
||||
|
||||
asm_opt.coverage = get_coverage(sources, coverage_cut, n_read);
|
||||
@@ -27154,9 +27168,11 @@ long long bubble_dist, int read_graph, int write)
|
||||
&R_INF, output_file_name);
|
||||
}
|
||||
|
||||
debug_info_of_specfic_read("m64062_190803_042216/120916077/ccs", sources, reverse_sources, -1, "beg");
|
||||
debug_info_of_specfic_read("m64062_190806_063919/71436446/ccs", sources, reverse_sources, -1, "beg");
|
||||
debug_info_of_specfic_read("m64062_190803_042216/82117654/ccs", sources, reverse_sources, -1, "beg");
|
||||
|
||||
|
||||
///debug_info_of_specfic_read("m64062_190803_042216/15205346/ccs", sources, reverse_sources, -1, "beg");
|
||||
|
||||
if (!(asm_opt.flag & HA_F_BAN_ASSEMBLY))
|
||||
{
|
||||
try_rescue_overlaps(sources, reverse_sources, n_read, 4);
|
||||
|
||||
@@ -690,6 +690,7 @@ void init_Debug_reads(Debug_reads* x, const char* file)
|
||||
x->query_num++;
|
||||
}
|
||||
x->read_name = (char**)malloc(sizeof(char*)*x->query_num);
|
||||
x->candidate_count = (kvec_t_u64_warp*)malloc(sizeof(kvec_t_u64_warp)*x->query_num);
|
||||
fseek(x->fp, 0, SEEK_SET);
|
||||
|
||||
i = 0;
|
||||
@@ -699,6 +700,7 @@ void init_Debug_reads(Debug_reads* x, const char* file)
|
||||
x->read_name[i] = (char*)malloc(sizeof(char)*(nameLen+1));
|
||||
memcpy(x->read_name[i], Name_Buffer, sizeof(char)*nameLen);
|
||||
x->read_name[i][nameLen] = '\0';
|
||||
kv_init(x->candidate_count[i].a);
|
||||
i++;
|
||||
}
|
||||
|
||||
@@ -716,6 +718,7 @@ void destory_Debug_reads(Debug_reads* x)
|
||||
for (i = 0; i < x->query_num; i++)
|
||||
{
|
||||
free(x->read_name[i]);
|
||||
kv_destroy(x->candidate_count[i].a);
|
||||
}
|
||||
|
||||
free(x->read_name);
|
||||
|
||||
@@ -148,6 +148,7 @@ typedef struct
|
||||
{
|
||||
char** read_name;
|
||||
uint64_t query_num;
|
||||
kvec_t_u64_warp* candidate_count;
|
||||
FILE* fp;
|
||||
pthread_mutex_t OutputMutex;
|
||||
} Debug_reads;
|
||||
|
||||
+5
-5
@@ -58,7 +58,7 @@ int ha_ov_type(const overlap_region *r, uint32_t len)
|
||||
}
|
||||
|
||||
void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag,
|
||||
kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar)
|
||||
kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct)
|
||||
{
|
||||
uint32_t i, rlen;
|
||||
uint64_t k, l;
|
||||
@@ -73,7 +73,7 @@ kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f
|
||||
rlen = Get_READ_LENGTH(R_INF, rid); // read length
|
||||
|
||||
// get the list of anchors
|
||||
ha_sketch_query(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag);
|
||||
ha_sketch_query(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag, dbg_ct);
|
||||
// minimizer of queried read
|
||||
if (ab->mz.m > ab->old_mz_m) {
|
||||
ab->old_mz_m = ab->mz.m;
|
||||
@@ -205,14 +205,14 @@ void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf)
|
||||
|
||||
|
||||
void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres,
|
||||
int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar)
|
||||
int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct)
|
||||
{
|
||||
extern void *ha_flt_tab;
|
||||
extern ha_pt_t *ha_idx;
|
||||
extern void *ha_flt_tab_hp;
|
||||
extern ha_pt_t *ha_idx_hp;
|
||||
|
||||
ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar);
|
||||
ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar, dbg_ct);
|
||||
|
||||
if(ha_idx_hp)
|
||||
{
|
||||
@@ -242,7 +242,7 @@ int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp*
|
||||
overlap_list->length = k;
|
||||
|
||||
|
||||
ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab_hp, ha_idx_hp, f_cigar);
|
||||
ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab_hp, ha_idx_hp, f_cigar, dbg_ct);
|
||||
|
||||
if(overlap_list->length + overlap_list_hp->length > overlap_list->size)
|
||||
{
|
||||
|
||||
@@ -37,6 +37,7 @@ void *ha_flt_tab;
|
||||
ha_pt_t *ha_idx;
|
||||
void *ha_flt_tab_hp;
|
||||
ha_pt_t *ha_idx_hp;
|
||||
void *ha_ct_table;
|
||||
|
||||
/***************************
|
||||
* Yak specific parameters *
|
||||
@@ -132,8 +133,8 @@ typedef struct {
|
||||
ha_ct1_t *h;
|
||||
} ha_ct_t;
|
||||
|
||||
///for 0-th counting, k = 51, pre = 12, n_hash = 4, n_shift = 0
|
||||
///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 37
|
||||
///for 0-th counting, k = 51, pre = 12, n_hash = 4, n_shift = 37
|
||||
///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0
|
||||
static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
|
||||
{
|
||||
ha_ct_t *h;
|
||||
@@ -146,8 +147,7 @@ static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
|
||||
///it seems there is a large hash table h, consisting 4096 small hash tables
|
||||
for (i = 0; i < 1<<h->pre; ++i)
|
||||
h->h[i].h = yak_ct_init();
|
||||
///for 0-th counting, don't enter here
|
||||
///seems used for minimzer
|
||||
///for 0-th counting, enter here; used for bloom filter
|
||||
if (n_hash > 0 && n_shift > h->pre) {
|
||||
h->n_hash = n_hash, h->n_shift = n_shift;
|
||||
for (i = 0; i < 1<<h->pre; ++i)
|
||||
@@ -749,10 +749,10 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
|
||||
} else if (c0) {
|
||||
pl.ct = c0, pl.create_new = !!(flag&HAF_CREATE_NEW);
|
||||
assert(c0->k == opt->k && c0->pre == opt->pre);
|
||||
} else {///for 0-th counting and 1-th counting, go into here
|
||||
} else {///for ft-th counting and 1-th counting, go into here
|
||||
pl.create_new = 1; // alware create new elements if the count table is empty
|
||||
///for 0-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0
|
||||
///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 37
|
||||
///for 0-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 37
|
||||
///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0
|
||||
///building a large hash table consisting of 4096 small hash tables
|
||||
pl.ct = ha_ct_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift);
|
||||
}
|
||||
@@ -785,11 +785,10 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
|
||||
opt.k = asm_opt->k_mer_length;
|
||||
///always 0
|
||||
opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC);
|
||||
///for 0-th counting, shoud be 1
|
||||
///for 1-th counting, shoud be 51
|
||||
///for ft-counting, shoud be 1
|
||||
opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win;
|
||||
///for 0-th counting, shoud be 0
|
||||
///for 1-th counting, shoud be 37
|
||||
///for ft-counting, shoud be 37
|
||||
///for ha_pt_gen, shoud be 0
|
||||
opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift;
|
||||
opt.n_thread = asm_opt->thread_num;
|
||||
///asm_opt->num_reads is the number of fastq files
|
||||
@@ -839,6 +838,32 @@ void ha_ft_destroy(void *h)
|
||||
if (h) yak_ft_destroy((yak_ft_t*)h);
|
||||
}
|
||||
|
||||
|
||||
void debug_ct_index(void* q_ct_idx, void* r_ct_idx)
|
||||
{
|
||||
ha_ct_t* ct_idx = (ha_ct_t*)q_ct_idx;
|
||||
yak_ct_t *g = NULL;
|
||||
uint64_t i;
|
||||
khint_t k;
|
||||
for (i = 0; (int)i < 1<<ct_idx->pre; i++)
|
||||
{
|
||||
g = ct_idx->h[i].h;
|
||||
for (k = 0; k < kh_end(g); ++k)
|
||||
{
|
||||
if (kh_exist(g, k))
|
||||
{
|
||||
int c = kh_key(g, k) & YAK_MAX_COUNT;
|
||||
uint64_t hash = ((kh_key(g, k) >> ct_idx->pre)<<ct_idx->pre) | i;
|
||||
int q = query_ct_index(r_ct_idx, hash);
|
||||
if(q!=c)
|
||||
{
|
||||
fprintf(stderr, "ERROR:c: %d, q: %d\n", c, q);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*************************
|
||||
* High-level interfaces *
|
||||
*************************/
|
||||
@@ -851,6 +876,15 @@ 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|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs);
|
||||
if((asm_opt->flag & HA_F_VERBOSE_GFA))
|
||||
{
|
||||
write_ct_index((void*)h, asm_opt->output_file_name);
|
||||
// load_ct_index(&ha_ct_table, asm_opt->output_file_name);
|
||||
// debug_ct_index((void*)h, ha_ct_table);
|
||||
// debug_ct_index(ha_ct_table, (void*)h);
|
||||
// ha_ct_destroy((ha_ct_t *)ha_ct_table);
|
||||
}
|
||||
|
||||
if(!(ex_flag & HAF_SKIP_READ))
|
||||
{
|
||||
ha_ct_hist(h, cnt, asm_opt->thread_num);
|
||||
@@ -917,9 +951,85 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
|
||||
return pt;
|
||||
}
|
||||
|
||||
int query_ct_index(void* ct_idx, uint64_t hash)
|
||||
{
|
||||
ha_ct1_t *g = &(((ha_ct_t*)ct_idx)->h[hash & ((1ULL<<((ha_ct_t*)ct_idx)->pre) - 1)]);
|
||||
khint_t k;
|
||||
k = yak_ct_get(g->h, hash);
|
||||
if (k == kh_end(g->h)) return 0;
|
||||
return kh_key(g->h, k)&YAK_MAX_COUNT;
|
||||
}
|
||||
|
||||
|
||||
int write_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
|
||||
int write_ct_index(void *i_ct_idx, char* file_name)
|
||||
{
|
||||
char* gfa_name = (char*)malloc(strlen(file_name)+25);
|
||||
sprintf(gfa_name, "%s.ct_flt", file_name);
|
||||
FILE* fp = fopen(gfa_name, "w");
|
||||
if (!fp) {
|
||||
free(gfa_name);
|
||||
return 0;
|
||||
}
|
||||
ha_ct_t* ct_idx = (ha_ct_t*)i_ct_idx;
|
||||
int i;
|
||||
ha_ct1_t *g;
|
||||
fwrite(&ct_idx->k, sizeof(ct_idx->k), 1, fp);
|
||||
fwrite(&ct_idx->pre, sizeof(ct_idx->pre), 1, fp);
|
||||
fwrite(&ct_idx->n_hash, sizeof(ct_idx->n_hash), 1, fp);
|
||||
fwrite(&ct_idx->n_shift, sizeof(ct_idx->n_shift), 1, fp);
|
||||
fwrite(&ct_idx->tot, sizeof(ct_idx->tot), 1, fp);
|
||||
for (i = 0; i < 1<<ct_idx->pre; i++)
|
||||
{
|
||||
g = &(ct_idx->h[i]);
|
||||
yak_ct_save(g->h, fp);
|
||||
}
|
||||
|
||||
|
||||
fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
|
||||
free(gfa_name);
|
||||
fclose(fp);
|
||||
return 1;
|
||||
}
|
||||
|
||||
int load_ct_index(void **i_ct_idx, char* file_name)
|
||||
{
|
||||
char* gfa_name = (char*)malloc(strlen(file_name)+25);
|
||||
sprintf(gfa_name, "%s.ct_flt", file_name);
|
||||
FILE* fp = fopen(gfa_name, "r");
|
||||
if (!fp) {
|
||||
free(gfa_name);
|
||||
return 0;
|
||||
}
|
||||
ha_ct_t** ct_idx = (ha_ct_t**)i_ct_idx;
|
||||
double index_time = 0;
|
||||
int i;
|
||||
ha_ct_t *h = 0;
|
||||
ha_ct1_t *g;
|
||||
CALLOC(h, 1);
|
||||
|
||||
fread(&h->k, sizeof(h->k), 1, fp);
|
||||
fread(&h->pre, sizeof(h->pre), 1, fp);
|
||||
fread(&h->n_hash, sizeof(h->n_hash), 1, fp);
|
||||
fread(&h->n_shift, sizeof(h->n_shift), 1, fp);
|
||||
fread(&h->tot, sizeof(h->tot), 1, fp);
|
||||
CALLOC(h->h, 1<<h->pre);
|
||||
|
||||
|
||||
index_time = yak_realtime();
|
||||
for (i = 0; i < 1<<h->pre; ++i)
|
||||
{
|
||||
g = &(h->h[i]);
|
||||
yak_ct_load(&(g->h), fp);
|
||||
}
|
||||
|
||||
(*ct_idx) = h;
|
||||
fprintf(stderr, "[M::%s::%.3f] ==> Loaded count table\n", __func__, yak_realtime() - index_time);
|
||||
fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__);
|
||||
free(gfa_name);
|
||||
return 1;
|
||||
}
|
||||
|
||||
int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
|
||||
{
|
||||
char* gfa_name = (char*)malloc(strlen(file_name)+25);
|
||||
sprintf(gfa_name, "%s.pt_flt", file_name);
|
||||
@@ -970,7 +1080,7 @@ int write_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt
|
||||
return 1;
|
||||
}
|
||||
|
||||
int load_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
|
||||
int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
|
||||
{
|
||||
char* gfa_name = (char*)malloc(strlen(file_name)+25);
|
||||
sprintf(gfa_name, "%s.pt_flt", file_name);
|
||||
|
||||
@@ -29,6 +29,8 @@ extern void *ha_flt_tab;
|
||||
extern ha_pt_t *ha_idx;
|
||||
extern void *ha_flt_tab_hp;
|
||||
extern ha_pt_t *ha_idx_hp;
|
||||
extern void *ha_ct_table;
|
||||
|
||||
|
||||
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode);
|
||||
int ha_ft_isflt(const void *hh, uint64_t y);
|
||||
@@ -38,8 +40,11 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
|
||||
void ha_pt_destroy(ha_pt_t *h);
|
||||
const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);
|
||||
|
||||
int write_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
|
||||
int load_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
|
||||
int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
|
||||
int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
|
||||
int write_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);
|
||||
|
||||
ha_abuf_t *ha_abuf_init(void);
|
||||
void ha_abuf_destroy(ha_abuf_t *ab);
|
||||
@@ -55,7 +60,7 @@ double yak_cpu_usage(void);
|
||||
void ha_triobin(const hifiasm_opt_t *opt);
|
||||
|
||||
void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf);
|
||||
void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag);
|
||||
void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct);
|
||||
int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het);
|
||||
|
||||
static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function
|
||||
|
||||
+5
-1
@@ -139,12 +139,15 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
|
||||
|
||||
|
||||
void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf,
|
||||
kvec_t_u8_warp* k_flag)
|
||||
kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct)
|
||||
{ ///in default, w = 51, k = 51, is_hpc = 1
|
||||
/**
|
||||
uint64_t x;
|
||||
uint64_t rid:28, pos:27, rev:1, span:8;
|
||||
**/
|
||||
extern void *ha_ct_table;
|
||||
if(dbg_ct != NULL) dbg_ct->a.n = 0;
|
||||
|
||||
static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 };
|
||||
uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0}, filtered;
|
||||
int i, j, l, buf_pos, min_pos, kmer_span = 0;
|
||||
@@ -201,6 +204,7 @@ kvec_t_u8_warp* k_flag)
|
||||
|
||||
filtered = 0;
|
||||
if(hf != 0) filtered = ha_ft_isflt(hf, y);
|
||||
if(dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));
|
||||
///if (hf == 0 || ha_ft_isflt(hf, y) == 0)
|
||||
if(filtered == 0)
|
||||
info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
|
||||
|
||||
Reference in New Issue
Block a user