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

save bfore multiple overlaps

Browse Source
This commit is contained in:
chhylp123
2020-08-23 08:46:33 -04:00
parent 90d3cb69b9
commit 03cf0d475c
15 changed files with 1120 additions and 131 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+193 -41
View File
@@ -11,10 +11,12 @@
#include "htab.h"
#include "kthread.h"
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);
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, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf);
void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list);
All_reads R_INF;
Debug_reads R_INF_FLAG;
void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int* new_length)
{
@@ -442,6 +444,7 @@ typedef struct {
UC_Read self_read, ovlp_read;
Candidates_list clist;
overlap_region_alloc olist;
overlap_region_alloc olist_hp;
ha_abuf_t *ab;
// error correction related buffers
int64_t num_read_base, num_correct_base, num_recorrect_base;
@@ -451,6 +454,9 @@ typedef struct {
Correct_dumy correct;
haplotype_evdience_alloc hap;
Round2_alignment round2;
kvec_t_u32_warp b_buf;
kvec_t_u64_warp r_buf;
kvec_t_u8_warp k_flag;
} ha_ovec_buf_t;
ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
@@ -462,6 +468,10 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
init_UC_Read(&b->ovlp_read);
init_Candidates_list(&b->clist);
init_overlap_region_alloc(&b->olist);
init_overlap_region_alloc(&b->olist_hp);
kv_init(b->b_buf.a);
kv_init(b->r_buf.a);
kv_init(b->k_flag.a);
b->ab = ha_abuf_init();
if (!b->is_final) {
init_Cigar_record(&b->cigar1);
@@ -480,7 +490,11 @@ void ha_ovec_destroy(ha_ovec_buf_t *b)
destory_UC_Read(&b->ovlp_read);
destory_Candidates_list(&b->clist);
destory_overlap_region_alloc(&b->olist);
destory_overlap_region_alloc(&b->olist_hp);
ha_abuf_destroy(b->ab);
kv_destroy(b->b_buf.a);
kv_destroy(b->r_buf.a);
kv_destroy(b->k_flag.a);
if (!b->is_final) {
destory_Cigar_record(&b->cigar1);
destory_Graph(&b->POA_Graph);
@@ -510,6 +524,7 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b)
{
int64_t i, mem = 0, mem_clist, mem_olist;
mem_clist = b->clist.size * sizeof(k_mer_hit) + b->clist.chainDP.size * 7 * 4;
mem_olist = b->olist.size * sizeof(overlap_region);
for (i = 0; i < (int64_t)b->olist.size; ++i) {
const overlap_region *r = &b->olist.list[i];
@@ -517,6 +532,14 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b)
mem_olist += r->f_cigar.size * 8;
mem_olist += r->boundary_cigars.size * sizeof(window_list);
}
mem_olist += b->olist_hp.size * sizeof(overlap_region);
for (i = 0; i < (int64_t)b->olist_hp.size; ++i) {
const overlap_region *r = &b->olist_hp.list[i];
mem_olist += r->w_list_size * sizeof(window_list);
mem_olist += r->f_cigar.size * 8;
mem_olist += r->boundary_cigars.size * sizeof(window_list);
}
mem = ha_abuf_mem(b->ab) + mem_clist + mem_olist;
if (!b->is_final) {
mem += sizeof(Cigar_record) + b->cigar1.lost_base_size + b->cigar1.size * 4;
@@ -534,8 +557,8 @@ static void worker_ovec(void *data, long i, int tid)
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
int fully_cov, abnormal;
ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 1);
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), &(R_INF.paf[i]), &(R_INF.reverse_paf[i]));
clear_Cigar_record(&b->cigar1);
clear_Round2_alignment(&b->round2);
@@ -557,7 +580,13 @@ static void worker_ovec(void *data, long i, int tid)
}
R_INF.paf[i].is_abnormal = abnormal;
if (b->save_ov) {
R_INF.trio_flag[i] = AMBIGU;
if(ha_idx_hp == NULL)
{
R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &b->b_buf, &b->r_buf, &(b->k_flag), RESEED_HP_RATE, NULL);
}
if (R_INF.trio_flag[i] != AMBIGU || b->save_ov) {
int is_rev = (asm_opt.number_of_round % 2 == 0);
push_overlaps(&(R_INF.paf[i]), &b->olist, 1, &R_INF, is_rev);
push_overlaps(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, is_rev);
@@ -568,25 +597,97 @@ static void worker_ovec(void *data, long i, int tid)
static void worker_ovec_related_reads(void *data, long i, int tid)
{
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
int required_read_name_length = strlen(asm_opt.required_read_name);
uint64_t k;
if (required_read_name_length == (int)Get_NAME_LENGTH((R_INF),i)
&&
memcmp(asm_opt.required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
{
ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 1);
fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i));
recover_UC_Read(&b->self_read, &R_INF, i);
fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq);
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)
{
break;
}
}
for (k = 0; k < b->olist.length; k++) {
fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
recover_UC_Read(&b->self_read, &R_INF, b->olist.list[k].y_id);
fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq);
if(k < R_INF_FLAG.query_num)
{
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), &(R_INF.paf[i]), &(R_INF.reverse_paf[i]));
clear_Cigar_record(&b->cigar1);
clear_Round2_alignment(&b->round2);
correct_overlap(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
&b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal);
b->num_read_base += b->self_read.length;
b->num_correct_base += b->correct.corrected_base;
b->num_recorrect_base += b->round2.dumy.corrected_base;
push_cigar(R_INF.cigars, i, &b->cigar1);
push_cigar(R_INF.second_round_cigar, i, &b->round2.cigar);
R_INF.paf[i].is_fully_corrected = 0;
if (fully_cov) {
if (get_cigar_errors(&b->cigar1) == 0 && get_cigar_errors(&b->round2.cigar) == 0)
R_INF.paf[i].is_fully_corrected = 1;
}
R_INF.paf[i].is_abnormal = abnormal;
pthread_mutex_lock(&R_INF_FLAG.OutputMutex);
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, "***************************forward ovlp***************************\n");
for (k = 0; k < b->olist.length; k++)
{
if(b->olist.list[k].is_match != 1) continue;
fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e,
b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
}
}
fprintf(R_INF_FLAG.fp, "***************************reverse ovlp***************************\n");
for (k = 0; k < b->olist.length; k++)
{
if(b->olist.list[k].is_match != 2) continue;
fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e,
b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
}
fprintf(R_INF_FLAG.fp, "***************************unmatched ovlp***************************\n");
for (k = 0; k < b->olist.length; k++)
{
if(b->olist.list[k].is_match == 1) continue;
if(b->olist.list[k].is_match == 2) continue;
fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e,
b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
}
R_INF.trio_flag[i] = AMBIGU;
if(ha_idx_hp == NULL)
{
R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &b->b_buf, &b->r_buf, &(b->k_flag), RESEED_HP_RATE, R_INF_FLAG.fp);
}
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);
}
}
@@ -693,25 +794,75 @@ void Output_corrected_reads()
fclose(output_file);
}
void debug_print_pob_regions()
{
uint64_t i, total = 0;
for (i = 0; i < R_INF.total_reads; i++)
{
if(R_INF.trio_flag[i]!=AMBIGU)
{
total++;
fprintf(stderr, "(%lu) %.*s\n", i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
}
}
fprintf(stderr, "total hp reads: %lu, R_INF.total_reads: %lu\n", total, R_INF.total_reads);
exit(1);
}
void rescue_hp_reads(ha_ovec_buf_t **b)
{
int hom_cov, het_cov;
ha_flt_tab_hp = ha_idx_hp = NULL;
if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
ha_flt_tab_hp = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 1);
}
ha_idx_hp = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 1, &R_INF, &hom_cov, &het_cov);
if (asm_opt.required_read_name)
kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads);
else
kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads);
ha_ft_destroy(ha_flt_tab_hp); ha_flt_tab_hp = NULL;
ha_pt_destroy(ha_idx_hp); ha_idx_hp = NULL;
}
void ha_overlap_and_correct(int round)
{
int i, hom_cov, het_cov;
int i, hom_cov, het_cov, r_out = 0;
ha_ovec_buf_t **b;
ha_ecsave_buf_t *e;
ha_flt_tab_hp = ha_idx_hp = NULL;
if((ha_idx == NULL)&&(asm_opt.flag & HA_F_VERBOSE_GFA)&&(round == asm_opt.number_of_round - 1))
{
r_out = 1;
}
if(asm_opt.required_read_name) init_Debug_reads(&R_INF_FLAG, asm_opt.required_read_name); // for debugging only
// overlap and correct reads
CALLOC(b, asm_opt.thread_num);
for (i = 0; i < asm_opt.thread_num; ++i)
b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1));
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF, &hom_cov, &het_cov); // build the index
if(ha_idx) hom_cov = asm_opt.hom_cov;
if(ha_idx == NULL) ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
ha_opt_update_cov(&asm_opt, hom_cov);
if (asm_opt.required_read_name)
kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads);
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);
ha_pt_destroy(ha_idx);
ha_idx = 0;
ha_idx = NULL;
// collect statistics
for (i = 0; i < asm_opt.thread_num; ++i) {
@@ -723,8 +874,8 @@ void ha_overlap_and_correct(int round)
}
free(b);
if (asm_opt.required_read_name) exit(0); // for debugging only
if (asm_opt.required_read_name) destory_Debug_reads(&R_INF_FLAG), exit(0); // for debugging only
// save corrected reads to R_INF
CALLOC(e, asm_opt.thread_num);
for (i = 0; i < asm_opt.thread_num; ++i) {
@@ -740,6 +891,7 @@ void ha_overlap_and_correct(int round)
free(e[i].second_round_read);
}
free(e);
debug_print_pob_regions();
}
@@ -1138,15 +1290,8 @@ static void worker_ov_final(void *data, long i, int tid)
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.001, asm_opt.max_n_chain, 0);
///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 0);
/**
correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon,
&matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1,
&current_cigar, &hap, &second_round, 0, 0);
push_final_overlaps(&(R_INF.paf[i]), &overlap_list);
**/
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), &(R_INF.paf[i]), &(R_INF.reverse_paf[i]));
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);
@@ -1211,7 +1356,8 @@ 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_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, HIGH_HET_ERROR_RATE, asm_opt.max_n_chain, 1);
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), &(R_INF.paf[i]), &(R_INF.reverse_paf[i]));
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);
@@ -1410,7 +1556,7 @@ void hap_recalculate_peaks(char* output_file_name)
int hom_cov, het_cov;
// construct hash table for high occurrence k-mers
if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov);
ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0);
ha_opt_update_cov(&asm_opt, hom_cov);
}
free(R_INF.read_length);
@@ -1418,7 +1564,7 @@ void hap_recalculate_peaks(char* output_file_name)
load_All_reads(&R_INF, gfa_name);
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
asm_opt.hom_cov = hom_cov;
asm_opt.het_cov = het_cov;
ha_pt_destroy(ha_idx);
@@ -1435,10 +1581,12 @@ void ha_overlap_final(void)
{
int i, hom_cov, het_cov;
ha_ovec_buf_t **b;
ha_flt_tab_hp = ha_idx_hp = NULL;
CALLOC(b, asm_opt.thread_num);
for (i = 0; i < asm_opt.thread_num; ++i)
b[i] = ha_ovec_init(asm_opt.flag & HA_F_HIGH_HET, 1);///b[i] = ha_ovec_init(1, 1);
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
if(asm_opt.flag & HA_F_HIGH_HET)
{
kt_for(asm_opt.thread_num, worker_ov_final_high_het, b, R_INF.total_reads);
@@ -1475,14 +1623,18 @@ int ha_assemble(void)
if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2;
}
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);
// construct hash table for high occurrence k-mers
if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov);
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);
ha_opt_update_cov(&asm_opt, hom_cov);
}
// error correction
assert(asm_opt.number_of_round > 0);
for (r = 0; r < asm_opt.number_of_round; ++r) {
for (r = ha_idx?asm_opt.number_of_round-1:0; 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
ha_overlap_and_correct(r);
fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(),
Assembly.h
+5
View File
@@ -8,6 +8,11 @@
#define Get_Cigar_Type(RECORD) (RECORD&3)
#define Get_Cigar_Length(RECORD) (RECORD>>2)
#define RESEED_DP 4
#define RESEED_PEAK_RATE 0.15
#define RESEED_LEN 2000
#define RESEED_HP_RATE 0.9
int ha_assemble(void);
#endif
CommandLines.cpp
+1 -1
View File
@@ -324,7 +324,7 @@ int check_option(hifiasm_opt_t* asm_opt)
if(asm_opt->fn_bin_yak[1] != NULL && check_file(asm_opt->fn_bin_yak[1], "YAK2") == 0) return 0;
if(asm_opt->fn_bin_list[0] != NULL && check_file(asm_opt->fn_bin_list[0], "LIST1") == 0) return 0;
if(asm_opt->fn_bin_list[1] != NULL && check_file(asm_opt->fn_bin_list[1], "LIST2") == 0) return 0;
if(asm_opt->required_read_name != NULL && check_file(asm_opt->required_read_name, "b") == 0) return 0;
// fprintf(stderr, "input file num: %d\n", asm_opt->num_reads);
// fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
// fprintf(stderr, "number of threads: %d\n", asm_opt->thread_num);
CommandLines.h
+1 -1
View File
@@ -3,7 +3,7 @@
#include <pthread.h>
#define HA_VERSION "0.11-r302"
#define HA_VERSION "0.11-r303"
#define VERBOSE 0
Correct.cpp
+103 -1
View File
@@ -7,8 +7,10 @@
#include "Assembly.h"
#include "CommandLines.h"
#include "ksw2.h"
#include "ksort.h"
#define generic_key(x) (x)
KRADIX_SORT_INIT(b32, uint32_t, generic_key, 4)
void clear_Round2_alignment(Round2_alignment* h)
@@ -7202,7 +7204,107 @@ void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R
}
void collect_no_cov_regions(overlap_region_alloc* overlap_list, All_reads* R_INF,
kvec_t_u32_warp* b, kvec_t_u64_warp* r, int min_dp, int min_len)
{
b->a.n = r->a.n = 0;
///if(overlap_list->length == 0) return;
long long i = 0, xLen = Get_READ_LENGTH((*R_INF), overlap_list->list[0].x_id);
uint32_t qs, qe;
uint64_t tmp;
int dp, old_dp, s_start = 0, s_end = 0;
///at least 1
if(min_len < 1) min_len = 1;
for (i = 0; i < (long long)overlap_list->length; i++)
{
if (overlap_list->list[i].is_match != 1 && overlap_list->list[i].is_match != 2) continue;
qs = overlap_list->list[i].x_pos_s;
qe = overlap_list->list[i].x_pos_e + 1;
kv_push(uint32_t, b->a, qs<<1);
kv_push(uint32_t, b->a, qe<<1|1);
}
///we can identify the qs and qe by the 0-th bit
radix_sort_b32(b->a.a, b->a.a + b->a.n);
for (i = 0, dp = 0; i < (long long)b->a.n; ++i)
{
old_dp = dp;
//if a[j] is qe
if (b->a.a[i]&1) --dp;
else ++dp;
/**
min_dp is the coverage drop threshold
there are two cases:
1. old_dp = dp + 1 (b.a[j] is qe); 2. old_dp = dp - 1 (b.a[j] is qs);
**/
if (old_dp < min_dp && dp >= min_dp) ///old_dp < dp, b.a[j] is qs
{
///case 2, a[j] is qs
s_end = b->a.a[i]>>1;
///at least 1
if(s_end-s_start >= min_len)
{
tmp = s_start; tmp = tmp << 32; tmp = tmp | (uint64_t)(s_end-1);
kv_push(uint64_t, r->a, tmp);
}
}
else if (old_dp >= min_dp && dp < min_dp) ///old_dp > min_dp, b.a[j] is qe
{
s_start = b->a.a[i]>>1;
}
}
if(s_start < xLen && xLen-s_start >= min_len)
{
s_end = xLen;
tmp = s_start; tmp = tmp << 32; tmp = tmp | (uint64_t)(s_end-1);
kv_push(uint64_t, r->a, tmp);
}
}
int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u32_warp* b, kvec_t_u64_warp* r, kvec_t_u8_warp* k_flag, float hp_rate, FILE* fp)
{
int i, k, qs, qe, ava_k_mer = 0, hp_k_mer = 0;
int min_dp = RESEED_DP;
if(asm_opt.hom_cov > 0) min_dp = asm_opt.hom_cov * RESEED_PEAK_RATE;
if(min_dp > RESEED_DP) min_dp = RESEED_DP;
collect_no_cov_regions(olist, R_INF, b, r, min_dp, RESEED_LEN);
for (i = 0; i < (int)r->a.n; i++)
{
///[qs, qe]
qs = r->a.a[i]>>32;
qe = (r->a.a[i]<<32)>>32;
for (k = qs; k <= qe; k++)
{
if(k_flag->a.a[k] > 1) ava_k_mer++;
if(k_flag->a.a[k] > 2) hp_k_mer++;
}
if(fp) fprintf(fp, "qs: %d, qe: %d, ava_k_mer: %d, hp_k_mer: %d\n", qs, qe, ava_k_mer, hp_k_mer);
}
if(fp) fprintf(fp, "ava_k_mer: %d, hp_k_mer: %d, hp_rate: %f\n", ava_k_mer, hp_k_mer, hp_rate);
if(fp)
{
for (k = 0; k < (int)k_flag->a.n; k++)
{
if(k_flag->a.a[k] > 0) fprintf(fp, "(%d) %u\n", k, k_flag->a.a[k]);
}
}
if(hp_k_mer > ava_k_mer*hp_rate) return 1; ///must use '>' instead of '>='
r->a.n = 0;
return 0;
}
void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read,
Correct.h
+1
View File
@@ -1173,6 +1173,7 @@ void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_I
UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read);
long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* overlap_read, uint8_t* x_num,
uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen);
int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u32_warp* b, kvec_t_u64_warp* r, kvec_t_u8_warp* k_flag, float hp_rate, FILE* fp);
#define FORWARD_KSW 0
#define BACKWARD_KSW 1
Hash_Table.cpp
+21
View File
@@ -391,6 +391,16 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
resize_Chain_Data(dp, a_n);
ret = ha_chain_check(a, a_n, dp, min_score, band_width_threshold);
/***************************************debug**************************************/
if(a_n > 0 && Get_NAME_LENGTH((R_INF),a[0].readID)==strlen("m64062_190803_042216/128778853/ccs"))
{
if (memcmp("m64062_190803_042216/128778853/ccs", Get_NAME((R_INF), a[0].readID),
Get_NAME_LENGTH((R_INF), a[0].readID)) == 0)
{
fprintf(stderr, "ret: %d\n", ret);
}
}
/***************************************debug**************************************/
if (ret > 0) {
a_n = ret;
goto skip_dp;
@@ -500,6 +510,17 @@ skip_dp:
}
/***************************************debug**************************************/
if(a_n > 0 && Get_NAME_LENGTH((R_INF),a[0].readID)==strlen("m64062_190803_042216/128778853/ccs"))
{
if (memcmp("m64062_190803_042216/128778853/ccs", Get_NAME((R_INF), a[0].readID),
Get_NAME_LENGTH((R_INF), a[0].readID)) == 0)
{
fprintf(stderr, "max_i: %lld\n", max_i);
}
}
/***************************************debug**************************************/
clear_fake_cigar(&(result->f_cigar));
///note a has been sorted by offset, that means has been sorted by query offset
i = max_i;
Overlaps.cpp
+134 -50
View File
@@ -2107,6 +2107,37 @@ 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)
{
fprintf(stderr, "\n\n\n");
uint32_t v, n_vtx = g->n_seq * 2, queryLen = strlen(name), flag = 0;
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);
if(g->seq[v>>1].del)
{
fprintf(stderr, "read has already been deleted.\n");
return;
}
asg_arc_t *av = asg_arc_a(g, v);
uint32_t i, nv = asg_arc_n(g, v);
for (i = 0; i < nv; ++i)
{
fprintf(stderr, "target: %.*s, el: %u, strong: %u, ol: %u, del: %u\n",
(int)Get_NAME_LENGTH(R_INF, (av[i].v>>1)),
Get_NAME(R_INF, (av[i].v>>1)),
av[i].el, av[i].strong, av[i].ol, av[i].del);
}
flag = 1;
}
}
}
asg_t *ma_sg_gen(const ma_hit_t_alloc* sources, long long n_read, const ma_sub_t *coverage_cut,
int max_hang, int min_ovlp)
@@ -8206,13 +8237,96 @@ int asg_arc_del_short_diploi_by_suspect_edge(asg_t *g, int max_ext)
return n_cut;
}
int if_potential_false_node(ma_hit_t_alloc* paf, uint64_t rLen, ma_sub_t* max_left, ma_sub_t* max_right, int if_set)
{
max_left->s = max_right->s = rLen;
max_left->e = max_right->e = 0;
long long j;
uint32_t qs, qe;
for (j = 0; j < paf->length; j++)
{
if(paf->buffer[j].del) continue;
if(paf->buffer[j].el != 1) continue;
qs = Get_qs(paf->buffer[j]);
qe = Get_qe(paf->buffer[j]);
///overlaps from left side
if(qs == 0)
{
if(qs < max_left->s) max_left->s = qs;
if(qe > max_left->e) max_left->e = qe;
}
///overlaps from right side
if(qe == rLen)
{
if(qs < max_right->s) max_right->s = qs;
if(qe > max_right->e) max_right->e = qe;
}
///note: if (qs == 0 && qe == rLen)
///this overlap would be added to both b_left and b_right
///that is what we want
}
if (max_left->e > max_right->s) return 0;
long long new_left_e, new_right_s;
new_left_e = max_left->e;
new_right_s = max_right->s;
for (j = 0; j < paf->length; j++)
{
if(paf->buffer[j].del) continue;
if(paf->buffer[j].el != 1) continue;
qs = Get_qs(paf->buffer[j]);
qe = Get_qe(paf->buffer[j]);
///check contained overlaps
if(qs != 0 && qe != rLen)
{
///[qs, qe), [max_left.s, max_left.e)
if(qs < max_left->e && qe > max_left->e)
{
///if(qe > max_left->e) max_left->e = qe;
if(qe > max_left->e && qe > new_left_e) new_left_e = qe;
}
///[qs, qe), [max_right.s, max_right.e)
if(qs < max_right->s && qe > max_right->s)
{
///if(qs < max_right->s) max_right->s = qs;
if(qs < max_right->s && qs < new_right_s) new_right_s = qs;
}
if(if_set)
{
max_left->e = new_left_e;
max_right->s = new_right_s;
}
}
}
max_left->e = new_left_e;
max_right->s = new_right_s;
if (max_left->e > max_right->s) return 0;
return 1;
}
//reomve edge between two chromesomes
//this node must be a single read
int asg_arc_del_false_node(asg_t *g, int max_ext)
int asg_arc_del_false_node(asg_t *g, ma_hit_t_alloc* sources, int max_ext)
{
double startTime = Get_T();
kvec_t(uint64_t) b;
memset(&b, 0, sizeof(b));
ma_sub_t max_left, max_right;
uint32_t v, n_vtx = g->n_seq * 2;
long long n_cut = 0;
@@ -8242,6 +8356,11 @@ int asg_arc_del_false_node(asg_t *g, int max_ext)
continue;
}
if(if_potential_false_node(&(sources[v>>1]), g->seq[v>>1].len, &max_left, &max_right, 1)==0)
{
continue;
}
asg_arc_t *av = asg_arc_a(g, v);
kv_push(uint64_t, b, (uint64_t)((uint64_t)av[0].ol << 32 | (av - g->arc)));
}
@@ -9582,35 +9701,6 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
void debug_info_of_specfic_node(const char* name, asg_t *g, char* command)
{
fprintf(stderr, "\n\n\n");
uint32_t v, n_vtx = g->n_seq * 2, queryLen = strlen(name);
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)
{
fprintf(stderr, "\nafter %s\n****************graph ref_read: %.*s, dir: %u****************\n",
command, (int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), v&1);
if(g->seq[v>>1].del)
{
fprintf(stderr, "read has already been deleted.\n");
return;
}
asg_arc_t *av = asg_arc_a(g, v);
uint32_t i, nv = asg_arc_n(g, v);
for (i = 0; i < nv; ++i)
{
fprintf(stderr, "target: %.*s, el: %u, strong: %u, ol: %u, del: %u\n",
(int)Get_NAME_LENGTH(R_INF, (av[i].v>>1)),
Get_NAME(R_INF, (av[i].v>>1)),
av[i].el, av[i].strong, av[i].ol, av[i].del);
}
return;
}
}
}
void debug_info_of_specfic_read(const char* name, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, int id, const char* command)
@@ -9636,8 +9726,8 @@ ma_hit_t_alloc* reverse_sources, int id, const char* command)
{
fprintf(stderr, "\n\n\nafter %s\n", command);
fprintf(stderr, "****************ma_hit_t (%lld)ref_read: %.*s****************\n",
i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
fprintf(stderr, "****************ma_hit_t (%lld)ref_read: %.*s, len: %lu****************\n",
i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), Get_READ_LENGTH(R_INF, i));
fprintf(stderr, "sources Len: %d, is_fully_corrected: %d\n",
@@ -12954,12 +13044,12 @@ float drop_ratio)
#define T_ROUND 2
asg_t *g = ug->g;
int round = T_ROUND;
redo:
asg_pop_bubble_primary_trio(ug, bubble_dist, (uint32_t)-1, DROP);
untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, DROP);
if(just_bubble_pop == 0)
{
cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex,
@@ -13013,7 +13103,7 @@ float drop_ratio)
{
if(round != T_ROUND)
{
unitig_arc_del_short_diploid_by_length(ug->g, drop_ratio);
unitig_arc_del_short_diploid_by_length(ug->g, drop_ratio);
}
round--;
goto redo;
@@ -22764,7 +22854,7 @@ kvec_asg_arc_t_warp* new_rtg_edges)
delete_useless_nodes(ug);
renew_utg(ug, read_g, new_rtg_edges);
}
if (!(asm_opt.flag & HA_F_BAN_POST_JOIN))
{
rescue_missing_overlaps_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang,
@@ -26792,9 +26882,6 @@ int max_hang, int min_ovlp)
}
}
void clean_graph(
int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
long long n_read, uint64_t* readLen, long long mini_overlap_length,
@@ -26815,7 +26902,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
///normalize_ma_hit_t_single_side(sources, n_read);
normalize_ma_hit_t_single_side_advance(sources, n_read);
normalize_ma_hit_t_single_side_advance(reverse_sources, n_read);
if (ha_opt_triobin(&asm_opt))
{
@@ -26840,11 +26926,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);
///debug_info_of_specfic_read("m54329U_190827_173812/166332272/ccs", sources, reverse_sources, -1, "clean");
sg = ma_sg_gen(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
asg_arc_del_trans(sg, gap_fuzz);
asm_opt.coverage = get_coverage(sources, coverage_cut, n_read);
if(VERBOSE >= 1)
@@ -26856,9 +26941,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
}
asg_cut_tip(sg, asm_opt.max_short_tip);
///debug_info_of_specfic_node("m64062_190803_042216/15205346/ccs", sg, "inner_1");
///drop_inexact_edegs_at_bubbles(sg, bubble_dist);
///debug_info_of_specfic_node("m54329U_190827_173812/166332272/ccs", sg, "beg");
if(clean_round > 0)
{
@@ -26896,7 +26980,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
asg_arc_identify_simple_bubbles_multi(sg, 1);
//reomve edge between two chromesomes
//this node must be a single read
asg_arc_del_false_node(sg, asm_opt.max_short_tip);
asg_arc_del_false_node(sg, sources, asm_opt.max_short_tip);
asg_cut_tip(sg, asm_opt.max_short_tip);
/****************************may have bugs********************************/
/****************************may have bugs********************************/
@@ -27044,10 +27128,7 @@ long long bubble_dist, int read_graph, int write)
init_R_to_U(&ruIndex, n_read);
asg_t *sg = NULL;
ma_sub_t* coverage_cut = NULL;
// debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", sources, reverse_sources, -1, "clean");
///actually min_thres = asm_opt.max_short_tip + 1 there are asm_opt.max_short_tip reads
min_thres = asm_opt.max_short_tip + 1;
@@ -27073,6 +27154,9 @@ long long bubble_dist, int read_graph, int write)
&R_INF, output_file_name);
}
///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);
Process_Read.cpp
+54
View File
@@ -46,6 +46,7 @@ void destory_All_reads(All_reads* r)
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->pb_regions) kv_destroy(r->pb_regions[i].a);
}
free(r->paf);
free(r->reverse_paf);
@@ -55,6 +56,7 @@ void destory_All_reads(All_reads* r)
free(r->name_index);
free(r->read_length);
free(r->trio_flag);
///if (r->pb_regions) free(r->pb_regions);
}
void write_All_reads(All_reads* r, char* read_file_name)
@@ -192,6 +194,7 @@ int load_All_reads(All_reads* r, char* read_file_name)
r->second_round_cigar[i].lost_base_length = r->cigars[i].lost_base_length = 0;
r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL;
}
///r->pb_regions = NULL;
free(index_name);
fclose(fp);
@@ -263,6 +266,8 @@ void malloc_All_reads(All_reads* r)
r->second_round_cigar = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
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++)
{
r->second_round_cigar[i].size = r->cigars[i].size = 0;
@@ -274,6 +279,7 @@ void malloc_All_reads(All_reads* r)
r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL;
init_ma_hit_t_alloc(&(r->paf[i]));
init_ma_hit_t_alloc(&(r->reverse_paf[i]));
///kv_init(r->pb_regions[i].a);
}
r->name = (char*)malloc(sizeof(char)*r->total_name_length);
@@ -667,3 +673,51 @@ void reverse_complement(char* pattern, uint64_t length)
pattern[end] = RC_CHAR(pattern[end]);
}
}
void init_Debug_reads(Debug_reads* x, const char* file)
{
int nameLen, i, bufLen = 1000;
if((uint64_t)(bufLen) < strlen(file) + 50) bufLen = strlen(file) + 50;
char* Name_Buffer = (char*)malloc(sizeof(char)*bufLen);
fprintf(stderr, "Queried debugging reads at: %s\n", file);
x->fp = fopen(file,"r");
x->query_num = 0;
while(fgets(Name_Buffer, bufLen, x->fp))
{
x->query_num++;
}
x->read_name = (char**)malloc(sizeof(char*)*x->query_num);
fseek(x->fp, 0, SEEK_SET);
i = 0;
while(fgets(Name_Buffer, bufLen, x->fp))
{
nameLen = strlen(Name_Buffer) - 1;
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';
i++;
}
fclose(x->fp);
sprintf(Name_Buffer, "%s.debug.stdout", file);
x->fp = fopen(Name_Buffer,"w");
fprintf(stderr, "Print debugging information to: %s\n", Name_Buffer);
free(Name_Buffer);
}
void destory_Debug_reads(Debug_reads* x)
{
uint64_t i;
for (i = 0; i < x->query_num; i++)
{
free(x->read_name[i]);
}
free(x->read_name);
fclose(x->fp);
}
Process_Read.h
+13 -1
View File
@@ -93,6 +93,7 @@ typedef struct
uint32_t new_length;
} Compressed_Cigar_record;
#define AMBIGU 0
#define FATHER 1
#define MOTHER 2
@@ -128,6 +129,8 @@ typedef struct
ma_hit_t_alloc* paf;
ma_hit_t_alloc* reverse_paf;
///kvec_t_u64_warp* pb_regions;
} All_reads;
extern All_reads R_INF;
@@ -140,6 +143,14 @@ typedef struct
long long RID;
} UC_Read;
typedef struct
{
char** read_name;
uint64_t query_num;
FILE* fp;
pthread_mutex_t OutputMutex;
} Debug_reads;
void init_All_reads(All_reads* r);
void malloc_All_reads(All_reads* r);
void ha_insert_read_len(All_reads *r, int read_len, int name_len);
@@ -154,5 +165,6 @@ void write_All_reads(All_reads* r, char* read_file_name);
int load_All_reads(All_reads* r, char* read_file_name);
void destory_All_reads(All_reads* r);
int destory_read_bin(All_reads* r);
void init_Debug_reads(Debug_reads* x, const char* file);
void destory_Debug_reads(Debug_reads* x);
#endif
anchor.cpp
+141 -7
View File
@@ -28,8 +28,8 @@ typedef struct {
} seed1_t;
struct ha_abuf_s {
uint64_t n_a, m_a;
uint32_t old_mz_m;
uint64_t n_a, m_a;///number of anchors (seed positions)
uint32_t old_mz_m;///number of seeds
ha_mz1_v mz;
seed1_t *seed;
anchor1_t *a;
@@ -57,10 +57,9 @@ static int ha_ov_type(const overlap_region *r, uint32_t len)
else return r->x_pos_s == 0? 0 : 1;
}
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)
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,
void *ha_flt_tab, ha_pt_t *ha_idx)
{
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
uint32_t i, rlen;
uint64_t k, l;
double low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
@@ -74,7 +73,8 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
rlen = Get_READ_LENGTH(R_INF, rid); // read length
// get the list of anchors
ha_sketch(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);
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);
// minimizer of queried read
if (ab->mz.m > ab->old_mz_m) {
ab->old_mz_m = ab->mz.m;
REALLOC(ab->seed, ab->old_mz_m);
@@ -93,6 +93,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
}
for (i = 0, k = 0; i < ab->mz.n; ++i) {
int j;
///z is one of the minimizer
ha_mz1_t *z = &ab->mz.a[i];
seed1_t *s = &ab->seed[i];
for (j = 0; j < s->n; ++j) {
@@ -116,6 +117,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
}
}
// copy over to _cl_
if (ab->m_a >= (uint64_t)cl->size) {
cl->size = ab->m_a;
@@ -128,6 +130,18 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
p->offset = ab->a[k].other_off;
p->self_offset = ab->a[k].self_off;
p->good = ab->a[k].good;
/***************************************debug**************************************/
if(Get_NAME_LENGTH((R_INF),p->readID)==strlen("m64062_190803_042216/128778853/ccs"))
{
if (memcmp("m64062_190803_042216/128778853/ccs", Get_NAME((R_INF), p->readID),
Get_NAME_LENGTH((R_INF), p->readID)) == 0)
{
fprintf(stderr, "(%lu) readID: %u, strand: %u, offset: %u, self_offset: %u\n",
k, p->readID, p->strand, p->offset, p->self_offset);
}
}
/***************************************debug**************************************/
}
cl->length = ab->n_a;
@@ -172,10 +186,130 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
}
}
ks_introsort_or_xs(overlap_list->length, overlap_list->list);
///ks_introsort_or_xs(overlap_list->length, overlap_list->list);
}
void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf)
{
uint64_t j = 0, inner_j = 0;
while (j < overlap_list->length && inner_j < paf->length)
{
if(overlap_list->list[j].y_id < paf->buffer[inner_j].tn)
{
j++;
}
else if(overlap_list->list[j].y_id > paf->buffer[inner_j].tn)
{
inner_j++;
}
else
{
if(overlap_list->list[j].y_pos_strand == paf->buffer[inner_j].rev)
{
overlap_list->list[j].is_match = 1;
}
j++;
inner_j++;
}
}
}
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, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf)
{
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, ha_flt_tab, ha_idx);
if(ha_idx_hp)
{
uint32_t i, k, y_id, overlapLen, max_i;
int shared_seed;
overlap_region t;
overlap_region_sort_y_id(overlap_list->list, overlap_list->length);
ma_hit_sort_tn(paf->buffer, paf->length);
ma_hit_sort_tn(rev_paf->buffer, rev_paf->length);
lable_matched_ovlp(overlap_list, paf);
lable_matched_ovlp(overlap_list, rev_paf);
for (i = 0, k = 0; i < overlap_list->length; ++i)
{
if(overlap_list->list[i].is_match == 1)
{
if(k != i)
{
t = overlap_list->list[k];
overlap_list->list[k] = overlap_list->list[i];
overlap_list->list[i] = t;
overlap_list->list[k].is_match = 0;
}
k++;
}
}
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, ha_flt_tab_hp, ha_idx_hp);
if(overlap_list->length + overlap_list_hp->length > overlap_list->size)
{
overlap_list->list = (overlap_region*)realloc(overlap_list->list,
sizeof(overlap_region)*(overlap_list->length + overlap_list_hp->length));
memset(overlap_list->list + overlap_list->size, 0, sizeof(overlap_region)*
(overlap_list->length + overlap_list_hp->length - overlap_list->size));
overlap_list->size = overlap_list->length + overlap_list_hp->length;
}
for (i = 0, k = overlap_list->length; i < overlap_list_hp->length; i++, k++)
{
t = overlap_list->list[k];
overlap_list->list[k] = overlap_list_hp->list[i];
overlap_list_hp->list[i] = t;
}
overlap_list->length = k;
overlap_region_sort_y_id(overlap_list->list, overlap_list->length);
i = k = 0;
while (i < overlap_list->length)
{
y_id = overlap_list->list[i].y_id;
shared_seed = overlap_list->list[i].shared_seed;
overlapLen = overlap_list->list[i].overlapLen;
max_i = i;
i++;
while (i < overlap_list->length && overlap_list->list[i].y_id == y_id)
{
if((overlap_list->list[i].shared_seed > shared_seed) ||
((overlap_list->list[i].shared_seed == shared_seed) && (overlap_list->list[i].overlapLen <= overlapLen)))
{
y_id = overlap_list->list[i].y_id;
shared_seed = overlap_list->list[i].shared_seed;
overlapLen = overlap_list->list[i].overlapLen;
max_i = i;
}
i++;
}
if(k != max_i)
{
t = overlap_list->list[k];
overlap_list->list[k] = overlap_list->list[max_i];
overlap_list->list[max_i] = t;
}
k++;
}
overlap_list->length = k;
}
ks_introsort_or_xs(overlap_list->length, overlap_list->list);
}
void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list)
{
htab.cpp
+262 -24
View File
@@ -35,6 +35,8 @@ const unsigned char seq_nt4_table[256] = { // translate ACGT to 0123
void *ha_flt_tab;
ha_pt_t *ha_idx;
void *ha_flt_tab_hp;
ha_pt_t *ha_idx_hp;
/***************************
* Yak specific parameters *
@@ -71,7 +73,7 @@ typedef struct {
int n_shift, n_hashes;
uint8_t *b;
} yak_bf_t;
///in most cases, n_shift = 25, n_hashes = 4
yak_bf_t *yak_bf_init(int n_shift, int n_hashes)
{
yak_bf_t *b;
@@ -126,10 +128,12 @@ typedef struct {
typedef struct {
int k, pre, n_hash, n_shift;
uint64_t tot;
uint64_t tot; ///number of distinct k-mers
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
static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
{
ha_ct_t *h;
@@ -138,12 +142,16 @@ static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
CALLOC(h, 1);
h->k = k, h->pre = pre;
CALLOC(h->h, 1<<h->pre);
///i<h->pre = 4096
///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
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)
h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash);
h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash); ///h->n_shift = 37, h->pre = 12, h->n_hash = 4
}
return h;
}
@@ -173,15 +181,23 @@ static int ha_ct_insert_list(ha_ct_t *h, int create_new, int n, const uint64_t *
int j, mask = (1<<h->pre) - 1, n_ins = 0;
ha_ct1_t *g;
if (n == 0) return 0;
///corresponding small hash index
g = &h->h[a[0]&mask];
for (j = 0; j < n; ++j) {
int ins = 1, absent;
///x is a 64-bit word, h->pre=12
///all elements at a have the same low 12 bits
///so low 12 bits are not useful
uint64_t x = a[j] >> h->pre;
khint_t k;
if ((a[j]&mask) != (a[0]&mask)) continue;
if (create_new) {
///for 0-th counting, g->b = NULL
if (g->b)
ins = (yak_bf_insert(g->b, x) == h->n_hash);
///for 0-th counting, g->b = NULL
///x = the high 52 bits of a[j] + low 12 bits 0
///the low 12 bits are used for counting
if (ins) {
k = yak_ct_put(g->h, x << YAK_COUNTER_BITS | (g->b? 1 : 0), &absent);
if (absent) ++n_ins;
@@ -219,6 +235,8 @@ static void worker_ct_hist(void *data, long i, int tid) // callback for kt_for()
++cnt[kh_key(g, k)&YAK_MAX_COUNT];
}
///YAK_N_COUNTS is also 4096
///used for calculating k-mer histogram
static void ha_ct_hist(const ha_ct_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread)
{
hist_aux_t a;
@@ -226,6 +244,7 @@ static void ha_ct_hist(const ha_ct_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread
a.h = h;
memset(cnt, 0, YAK_N_COUNTS * sizeof(uint64_t));
CALLOC(a.cnt, n_thread);
///start 4096 threads
kt_for(n_thread, worker_ct_hist, &a, 1<<h->pre);
for (i = 0; i < YAK_N_COUNTS; ++i) cnt[i] = 0;
for (j = 0; j < n_thread; ++j)
@@ -265,6 +284,7 @@ static void ha_ct_shrink(ha_ct_t *h, int min, int max, int n_thread)
int i;
shrink_aux_t a;
a.h = h, a.min = min, a.max = max;
///still start 4096 threads
kt_for(n_thread, worker_ct_shrink, &a, 1<<h->pre);
for (i = 0, h->tot = 0; i < 1<<h->pre; ++i)
h->tot += kh_size(h->h[i].h);
@@ -306,6 +326,7 @@ static void worker_pt_gen(void *data, long i, int tid) // callback for kt_for()
int absent;
khint_t l;
l = yak_pt_put(b->h, kh_key(g, k) >> a->ct->pre << YAK_COUNTER_BITS, &absent);
///this should be the start index of kh_key's corresponding pos at ha_idxpos_t* a
kh_val(b->h, l) = b->n;
b->n += kh_key(g, k) & YAK_MAX_COUNT;
}
@@ -412,8 +433,12 @@ typedef struct {
ha_mz1_t *b;
} ch_buf_t;
///p = 12
static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer
{
///assign k-mer to one of the 4096 bins
///using low 12 bits for assigning
///so all elements at b have the same low 12 bits
int pre = y & ((1<<p) - 1);
ch_buf_t *b = &buf[pre];
if (b->n == b->m) {
@@ -425,6 +450,7 @@ static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k
static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y)
{
///assign minimizer to one of 4096 bins by low 12 bits
int pre = y->x & ((1<<p) - 1);
ch_buf_t *b = &buf[pre];
if (b->n == b->m) {
@@ -434,13 +460,17 @@ static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y)
b->b[b->n++] = *y;
}
///buf is the read block, k is the k-mer length, p = 12, len is the read length, seq is the read
static void count_seq_buf(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf
{
int i, l;
uint64_t x[4], mask = (1ULL<<k) - 1, shift = k - 1;
for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) {
int c = seq_nt4_table[(uint8_t)seq[i]];
///c = 00, 01, 10, 11
if (c < 4) { // not an "N" base
///x[0] & x[1] are the forward k-mer
///x[2] & x[3] are the reverse complementary k-mer
x[0] = (x[0] << 1 | (c&1)) & mask;
x[1] = (x[1] << 1 | (c>>1)) & mask;
x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift;
@@ -483,12 +513,13 @@ KSEQ_INIT(gzFile, gzread)
#define HAF_RS_WRITE_SEQ 0x8
#define HAF_RS_READ 0x10
#define HAF_CREATE_NEW 0x20
#define HAF_SKIP_READ 0x40
typedef struct { // global data structure for kt_pipeline()
const yak_copt_t *opt;
const void *flt_tab;
int flag, create_new, is_store;
uint64_t n_seq;
uint64_t n_seq; ///number of total reads
kseq_t *ks;
UC_Read ucr;
ha_ct_t *ct;
@@ -499,7 +530,8 @@ typedef struct { // global data structure for kt_pipeline()
typedef struct { // data structure for each step in kt_pipeline()
pl_data_t *p;
uint64_t n_seq0;
uint64_t n_seq0; ///the start index of current buffer block at R_INF
///sum_len = total bases, nk = number of k-mers
int n_seq, m_seq, sum_len, nk;
int *len;
char **seq;
@@ -514,15 +546,17 @@ static void worker_for_insert(void *data, long i, int tid) // callback for kt_fo
ch_buf_t *b = &s->buf[i];
if (s->p->pt)
b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b);
else
else///for 0-th count, go into here
b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a);
}
static void worker_for_mz(void *data, long i, int tid)
{
st_data_t *s = (st_data_t*)data;
///get the corresponding minimzer vector of this read
ha_mz1_v *b = &s->mz_buf[tid];
s->mz_buf[tid].n = 0;
///s->p->opt->w = 51, s->p->opt->k
ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab);
s->mz[i].n = s->mz[i].m = b->n;
MALLOC(s->mz[i].a, b->n);
@@ -540,6 +574,11 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
s->n_seq0 = p->n_seq;
if (p->rs_in && (p->flag & HAF_RS_READ)) {
while (p->n_seq < p->rs_in->total_reads) {
if((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU)
{
++p->n_seq;
continue;
}
int l;
recover_UC_Read(&p->ucr, p->rs_in, p->n_seq);
l = p->ucr.length;
@@ -565,6 +604,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
exit(1);
}
if (p->rs_out) {
///for 0-th count, just insert read length to R_INF, instead of read
if (p->flag & HAF_RS_WRITE_LEN) {
assert(p->n_seq == p->rs_out->total_reads);
ha_insert_read_len(p->rs_out, l, p->ks->name.l);
@@ -578,6 +618,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);
}
}
///for 0-th count, insert both seq and length to local block
if (s->n_seq == s->m_seq) {
s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
REALLOC(s->len, s->m_seq);
@@ -589,6 +630,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
++p->n_seq;
s->sum_len += l;
s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
///p->opt->chunk_size is the block max size
if (s->sum_len >= p->opt->chunk_size)
break;
}
@@ -596,18 +638,24 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
if (s->sum_len == 0) free(s);
else return s;
} else if (step == 1) { // step 2: extract k-mers
///s is the block of reads
st_data_t *s = (st_data_t*)in;
///for 0-th counting, n_pre = 4096
int i, n_pre = 1<<p->opt->pre, m;
// allocate the k-mer buffer
CALLOC(s->buf, n_pre);
m = (int)(s->nk * 1.2 / n_pre) + 1;
//pre-allocate memory for each of 4096 buffer
for (i = 0; i < n_pre; ++i) {
s->buf[i].m = m;
///for 0-th counting, p->pt = NULL
if (p->pt) MALLOC(s->buf[i].b, m);
else MALLOC(s->buf[i].a, m);
}
// fill the buffer
///for 0-th counting, p->opt->w == 1
if (p->opt->w == 1) { // enumerate all k-mers
///scan all reads
for (i = 0; i < s->n_seq; ++i) {
if (p->opt->is_HPC)
count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
@@ -618,18 +666,21 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
} else { // minimizers only
uint32_t j;
// compute minimizers
// 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);
///calculate minimzers for each read, each read corresponds to one thread
kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq);
for (i = 0; i < p->opt->n_thread; ++i)
free(s->mz_buf[i].a);
free(s->mz_buf);
// insert minimizers
if (p->pt) {
if (p->pt) {///insert whole minimizer
for (i = 0; i < s->n_seq; ++i)
for (j = 0; j < s->mz[i].n; ++j)
pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]);
} else {
} else {///just insert the hash key of minimizer
for (i = 0; i < s->n_seq; ++i)
for (j = 0; j < s->mz[i].n; ++j)
ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
@@ -640,6 +691,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
}
free(s->mz);
}
///just clean seq
free(s->seq); free(s->len);
s->seq = 0, s->len = 0;
return s;
@@ -647,7 +699,9 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
st_data_t *s = (st_data_t*)in;
int i, n = 1<<p->opt->pre;
uint64_t n_ins = 0;
///for 0-th counting, p->pt = NULL
kt_for(p->opt->n_thread, worker_for_insert, s, n);
///n_ins is number of distinct k-mers
for (i = 0; i < n; ++i) {
n_ins += s->buf[i].n_ins;
if (p->pt) free(s->buf[i].b);
@@ -668,6 +722,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
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, int64_t *n_seq)
{
///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN|HAF_CREATE_NEW
int read_rs = (rs && (flag & HAF_RS_READ));
pl_data_t pl;
gzFile fp = 0;
@@ -676,23 +731,29 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
if (read_rs) {
pl.rs_in = rs;
init_UC_Read(&pl.ucr);
} else {
} else {///for 0-th counting, go into here
if ((fp = gzopen(fn, "r")) == 0) return 0;
pl.ks = kseq_init(fp);
}
///for 0-th counting, read all reads into pl.rs_out
if (rs && (flag & (HAF_RS_WRITE_LEN|HAF_RS_WRITE_SEQ)))
pl.rs_out = rs;
///for 0-th counting, flt_tab = NULL
///for 1-th counting, flt_tab = NULL
pl.flt_tab = flt_tab;
pl.opt = opt;
pl.flag = flag;
if (p0) {
if (p0) {///for 1-th counting, p0 = NULL
pl.pt = p0, pl.create_new = 0; // never create new elements in a position table
assert(p0->k == opt->k && p0->pre == opt->pre);
} else if (c0) {
pl.ct = c0, pl.create_new = !!(flag&HAF_CREATE_NEW);
assert(c0->k == opt->k && c0->pre == opt->pre);
} else {
} else {///for 0-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
///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);
}
kt_pipeline(3, worker_count, &pl, 3);
@@ -713,6 +774,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
yak_copt_t opt;
ha_ct_t *h = 0;
assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both
///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN
if (rs) {
if (flag & HAF_RS_WRITE_LEN)
init_All_reads(rs);
@@ -721,10 +783,16 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
}
yak_copt_init(&opt);
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
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
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
for (i = 0; i < asm_opt->num_reads; ++i)
h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, &n_seq);
if (h && opt.bf_shift > 0)
@@ -775,19 +843,24 @@ void ha_ft_destroy(void *h)
* High-level interfaces *
*************************/
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode)
{
yak_ft_t *flt_tab;
int64_t cnt[YAK_N_COUNTS];
int peak_hom, peak_het, cutoff;
int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0;
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, NULL, NULL, rs);
ha_ct_hist(h, cnt, asm_opt->thread_num);
peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
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);
cutoff = (int)(peak_hom * asm_opt->high_factor);
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs);
if(!(ex_flag & HAF_SKIP_READ))
{
ha_ct_hist(h, cnt, asm_opt->thread_num);
peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
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);
///in default, asm_opt->high_factor = 5.0
cutoff = (int)(peak_hom * asm_opt->high_factor);
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
}
ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num);
flt_tab = gen_hh(h);
ha_ct_destroy(h);
@@ -796,21 +869,23 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
return (void*)flt_tab;
}
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov)
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov)
{
int64_t cnt[YAK_N_COUNTS], tot_cnt;
int peak_hom, peak_het, i, extra_flag1, extra_flag2;
ha_ct_t *ct;
ha_pt_t *pt;
if (read_from_store) {
if (read_from_store) {///if reads have already been read
extra_flag1 = extra_flag2 = HAF_RS_READ;
} else if (rs->total_reads == 0) {
} else if (rs->total_reads == 0) {///if reads & length have not been scanned
extra_flag1 = HAF_RS_WRITE_LEN;
extra_flag2 = HAF_RS_WRITE_SEQ;
} else {
} else {///if length has been loaded but reads have not
extra_flag1 = HAF_RS_WRITE_SEQ;
extra_flag2 = HAF_RS_READ;
}
if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ;
ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
yak_realtime(), yak_cpu_usage(), (long)ct->tot);
@@ -820,12 +895,16 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
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);
///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);
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
if((extra_flag1 & HAF_SKIP_READ) && (extra_flag2 & HAF_SKIP_READ)) cutoff = YAK_MAX_COUNT - 1;
ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num);
for (i = 2, tot_cnt = 0; i <= cutoff; ++i) tot_cnt += cnt[i] * i;
} else {
///Note: here is just to remove minimizer appearing YAK_MAX_COUNT times
///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen
ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num);
for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i;
}
@@ -837,3 +916,162 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos);
return pt;
}
int write_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);
FILE* fp = fopen(gfa_name, "w");
if (!fp) {
free(gfa_name);
return 0;
}
yak_ft_t *ha_flt_tab = (yak_ft_t*)flt_tab;
if(ha_flt_tab)
{
fwrite("f", 1, 1, fp);
yak_ft_save(ha_flt_tab, fp);
}
if(ha_idx)
{
int i;
ha_pt1_t *g;
fwrite("h", 1, 1, fp);
fwrite(&ha_idx->k, sizeof(ha_idx->k), 1, fp);
fwrite(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp);
fwrite(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp);
fwrite(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp);
for (i = 0; i < 1<<ha_idx->pre; ++i)
{
g = &(ha_idx->h[i]);
yak_pt_save(g->h, fp);
fwrite(&g->n, sizeof(g->n), 1, fp);
fwrite(g->a, sizeof(ha_idxpos_t), g->n, fp);
}
}
fwrite(&opt->number_of_round, sizeof(opt->number_of_round), 1, fp);
fwrite(&opt->hom_cov, sizeof(opt->hom_cov), 1, fp);
fwrite(&opt->het_cov, sizeof(opt->het_cov), 1, fp);
fwrite(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp);
write_All_reads(r, gfa_name);
fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
free(gfa_name);
fclose(fp);
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)
{
char* gfa_name = (char*)malloc(strlen(file_name)+25);
sprintf(gfa_name, "%s.pt_flt", file_name);
FILE* fp = fopen(gfa_name, "r");
if (!fp) {
free(gfa_name);
return 0;
}
ha_pt_t *ha_idx = NULL;
char mode = 0;
int f_flag, absent, i;
double index_time, index_s_time, pos_time, pos_s_time;
f_flag += fread(&mode, 1, 1, fp);
if(mode == 'f')
{
index_time = yak_realtime();
yak_ft_load((yak_ft_t **)r_flt_tab, fp);
f_flag += fread(&mode, 1, 1, fp);
fprintf(stderr, "[M::%s::%.3f] ==> Loaded flt table\n", __func__, yak_realtime()-index_time);
}
///insert using multiple threads???
if(mode == 'h')
{
pos_time = index_time = 0;
CALLOC(ha_idx, 1);
ha_pt1_t *g;
f_flag += fread(&ha_idx->k, sizeof(ha_idx->k), 1, fp);
f_flag += fread(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp);
f_flag += fread(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp);
f_flag += fread(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp);
CALLOC(ha_idx->h, 1<<ha_idx->pre);
for (i = 0; i < 1<<ha_idx->pre; ++i)
{
index_s_time = yak_realtime();
g = &(ha_idx->h[i]);
yak_pt_load(&(g->h), fp);
index_time += yak_realtime() - index_s_time;
pos_s_time = yak_realtime();
f_flag += fread(&g->n, sizeof(g->n), 1, fp);
MALLOC(g->a, g->n);
f_flag += fread(g->a, sizeof(ha_idxpos_t), g->n, fp);
pos_time += yak_realtime() - pos_s_time;
}
(*r_ha_idx) = ha_idx;
fprintf(stderr, "[M::%s::%.3f(index)/%.3f(pos)] ==> Loaded pos table\n", __func__, index_time, pos_time);
}
if(mode != 'h' && mode != 'f')
{
free(gfa_name);
fclose(fp);
return 0;
}
f_flag += fread(&absent, sizeof(absent), 1, fp);
if(absent != opt->number_of_round)
{
fprintf(stderr, "ERROR: different number of rounds!\n");
exit(1);
}
f_flag += fread(&opt->hom_cov, sizeof(opt->hom_cov), 1, fp);
f_flag += fread(&opt->het_cov, sizeof(opt->het_cov), 1, fp);
f_flag += fread(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp);
fclose(fp);
if(!load_All_reads(r, gfa_name))
{
free(gfa_name);
return 0;
}
memset(r->trio_flag, AMBIGU, r->total_reads*sizeof(uint8_t));
r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
for (i = 0; i < (long long)r->total_reads; i++)
{
init_ma_hit_t_alloc(&(r->paf[i]));
init_ma_hit_t_alloc(&(r->reverse_paf[i]));
}
fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__);
free(gfa_name);
return 1;
}
htab.h
+12 -3
View File
@@ -6,7 +6,9 @@
#include "CommandLines.h"
typedef struct {
uint64_t x;
uint64_t x; ///x is the hash key
///rid is the read id, pos is the end pos of this minimizer, rev is the direction
///span is the length of this k-mer. For non-HPC k-mer, span may not be equal to k
uint64_t rid:28, pos:27, rev:1, span:8;
} ha_mz1_t;
@@ -25,15 +27,20 @@ typedef struct ha_abuf_s ha_abuf_t;
extern const unsigned char seq_nt4_table[256];
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
extern void *ha_flt_tab_hp;
extern ha_pt_t *ha_idx_hp;
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov);
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);
void ha_ft_destroy(void *h);
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov);
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov);
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);
ha_abuf_t *ha_abuf_init(void);
void ha_abuf_destroy(ha_abuf_t *ab);
uint64_t ha_abuf_mem(const ha_abuf_t *ab);
@@ -48,6 +55,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);
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
@@ -76,6 +84,7 @@ static inline uint64_t yak_hash64_64(uint64_t key)
static inline uint64_t yak_hash_long(uint64_t x[4])
{
///compare forward k-mer and reverse complementary strand
int j = x[1] < x[3]? 0 : 1;
return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]);
}
khashl.h
+32
View File
@@ -133,6 +133,29 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
h->count = 0; \
} \
}
#define __KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \
SCOPE khint_t prefix##_save(HType *h, FILE* fp) { \
if (!h) return 0; \
khint_t n_buckets = (h->keys? 1U<<h->bits : 0U); \
fwrite(&n_buckets, sizeof(n_buckets), 1, fp); \
fwrite(&h->bits, sizeof(h->bits), 1, fp); \
fwrite(&h->count, sizeof(h->count), 1, fp); \
fwrite(h->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
fwrite(h->keys, sizeof(khkey_t), n_buckets, fp); \
return 1; \
} \
SCOPE khint_t prefix##_load(HType **h, FILE* fp) { \
(*h) = prefix##_init(); \
khint_t n_buckets; \
fread(&n_buckets, sizeof(n_buckets), 1, fp); \
fread(&(*h)->bits, sizeof((*h)->bits), 1, fp); \
fread(&(*h)->count, sizeof((*h)->count), 1, fp); \
(*h)->used = (khint32_t*)kmalloc(__kh_fsize(n_buckets) * sizeof(khint32_t)); \
(*h)->keys = (khkey_t*)kmalloc(n_buckets * sizeof(khkey_t)); \
fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \
return 1; \
} \
#define __KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
SCOPE khint_t prefix##_getp(const HType *h, const khkey_t *key) { \
@@ -245,6 +268,7 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
#define KHASHL_INIT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
__KHASHL_TYPE(HType, khkey_t) \
__KHASHL_IMPL_BASIC(SCOPE, HType, prefix) \
__KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \
__KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
__KHASHL_IMPL_RESIZE(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
__KHASHL_IMPL_PUT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
@@ -264,6 +288,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
KHASHL_INIT(KH_LOCAL, HType, prefix##_s, HType##_s_bucket_t, prefix##_s_hash, prefix##_s_eq) \
SCOPE HType *prefix##_init(void) { return prefix##_s_init(); } \
SCOPE void prefix##_destroy(HType *h) { prefix##_s_destroy(h); } \
SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_s_save(h, fp); } \
SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_s_load(h, fp); } \
SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_s_resize(h, new_n_buckets); } \
SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_s_bucket_t t; t.key = key; return prefix##_s_getp(h, &t); } \
SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_s_del(h, k); } \
@@ -276,6 +302,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
KHASHL_INIT(KH_LOCAL, HType, prefix##_m, HType##_m_bucket_t, prefix##_m_hash, prefix##_m_eq) \
SCOPE HType *prefix##_init(void) { return prefix##_m_init(); } \
SCOPE void prefix##_destroy(HType *h) { prefix##_m_destroy(h); } \
SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_m_save(h, fp); } \
SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_m_load(h, fp); } \
SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_m_resize(h, new_n_buckets); } \
SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_m_bucket_t t; t.key = key; return prefix##_m_getp(h, &t); } \
SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_m_del(h, k); } \
@@ -287,6 +315,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
KHASHL_INIT(KH_LOCAL, HType, prefix##_cs, HType##_cs_bucket_t, __kh_cached_hash, prefix##_cs_eq) \
SCOPE HType *prefix##_init(void) { return prefix##_cs_init(); } \
SCOPE void prefix##_destroy(HType *h) { prefix##_cs_destroy(h); } \
SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_cs_save(h, fp); } \
SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_cs_load(h, fp); } \
SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cs_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cs_getp(h, &t); } \
SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cs_del(h, k); } \
SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cs_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cs_putp(h, &t, absent); }
@@ -297,6 +327,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
KHASHL_INIT(KH_LOCAL, HType, prefix##_cm, HType##_cm_bucket_t, __kh_cached_hash, prefix##_cm_eq) \
SCOPE HType *prefix##_init(void) { return prefix##_cm_init(); } \
SCOPE void prefix##_destroy(HType *h) { prefix##_cm_destroy(h); } \
SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_cm_save(h, fp); } \
SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_cm_load(h, fp); } \
SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cm_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cm_getp(h, &t); } \
SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cm_del(h, k); } \
SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cm_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cm_putp(h, &t, absent); }
sketch.cpp
+147 -2
View File
@@ -37,7 +37,11 @@ static inline int tq_shift(tiny_queue_t *q)
* @param p minimizers
*/
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)
{
{ ///in default, w = 51, k = 51, is_hpc = 1
/**
uint64_t x;
uint64_t rid:28, pos:27, rev:1, span:8;
**/
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};
int i, j, l, buf_pos, min_pos, kmer_span = 0;
@@ -45,8 +49,10 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
tiny_queue_t tq;
assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
///sizeof(ha_mz1_t) = 16
memset(buf, 0xff, w * 16);
memset(&tq, 0, sizeof(tiny_queue_t));
///len/w is the evaluated minimizer numbers
kv_resize(ha_mz1_t, *p, p->n + len/w);
for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
@@ -65,7 +71,11 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
tq_push(&tq, skip_len);
kmer_span += skip_len;
if (tq.count > k) kmer_span -= tq_shift(&tq);
} else kmer_span = l + 1 < k? l + 1 : k;
} else kmer_span = l + 1 < k? l + 1 : k;
///kmer_span should be used for HPC k-mer
///so for non-HPC k-mer, kmer_span should be k in any case?
///kmer_span is used to calculate anchor pos on reverse complementary strand
kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer
kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer
@@ -80,6 +90,12 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
}
} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
//for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i
//i is the real base iterator, while l is the HPC base iterator
//only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX)
//but even if l < k, infor is still stored into buf
buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
for (j = buf_pos + 1; j < w; ++j)
@@ -87,15 +103,26 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
for (j = 0; j < buf_pos; ++j)
if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
}
/**
* There are three cases:
* 1. info.x <= min.x, means info is a new minimizer
* 2. info.x > min.x, info is not a new minimizer
* (1) buf_pos != min_pos, do nothing
* (2) buf_pos == min_pos, means current minimizer has moved outside the window
* **/
///three cases: 1.
if (info.x <= min.x) { // a new minimum; then write the old min
if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
min = info, min_pos = buf_pos;
} else if (buf_pos == min_pos) { // old min has moved outside the window
if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
///buf_pos == min_pos, means current minimizer has moved outside the window
///so for now we need to find a new minimizer at the current window (w k-mers)
for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers
if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
for (j = 0; j <= buf_pos; ++j)
if (min.x >= buf[j].x) min = buf[j], min_pos = j;
if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
@@ -108,3 +135,121 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
if (min.x != UINT64_MAX)
kv_push(ha_mz1_t, *p, min);
}
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)
{ ///in default, w = 51, k = 51, is_hpc = 1
/**
uint64_t x;
uint64_t rid:28, pos:27, rev:1, span:8;
**/
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;
ha_mz1_t buf[256], min = dummy;
tiny_queue_t tq;
if(k_flag != NULL)
{
kv_resize(uint8_t, k_flag->a, (uint64_t)len);
k_flag->a.n = len;
memset(k_flag->a.a, 0, k_flag->a.n);
}
assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
///sizeof(ha_mz1_t) = 16
memset(buf, 0xff, w * 16);
memset(&tq, 0, sizeof(tiny_queue_t));
///len/w is the evaluated minimizer numbers
kv_resize(ha_mz1_t, *p, p->n + len/w);
for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
int c = seq_nt4_table[(uint8_t)str[i]];
ha_mz1_t info = dummy;
if (c < 4) { // not an ambiguous base
int z;
if (is_hpc) {
int skip_len = 1;
if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {
for (skip_len = 2; i + skip_len < len; ++skip_len)
if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)
break;
i += skip_len - 1; // put $i at the end of the current homopolymer run
}
tq_push(&tq, skip_len);
kmer_span += skip_len;
if (tq.count > k) kmer_span -= tq_shift(&tq);
} else kmer_span = l + 1 < k? l + 1 : k;
///kmer_span should be used for HPC k-mer
///so for non-HPC k-mer, kmer_span should be k in any case?
///kmer_span is used to calculate anchor pos on reverse complementary strand
if(k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC
kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer
kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer
kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1;
if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand
z = kmer[1] < kmer[3]? 0 : 1; // strand
++l;
if (l >= k && kmer_span < 256) {
uint64_t y;
y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
filtered = 0;
if(hf != 0) filtered = ha_ft_isflt(hf, y);
///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;
if(k_flag != NULL) k_flag->a.a[i]++;
if(k_flag != NULL && filtered == 1) k_flag->a.a[i]++;
}
} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
//for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i
//i is the real base iterator, while l is the HPC base iterator
//only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX)
//but even if l < k, infor is still stored into buf
buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
for (j = buf_pos + 1; j < w; ++j)
if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
for (j = 0; j < buf_pos; ++j)
if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
}
/**
* There are three cases:
* 1. info.x <= min.x, means info is a new minimizer
* 2. info.x > min.x, info is not a new minimizer
* (1) buf_pos != min_pos, do nothing
* (2) buf_pos == min_pos, means current minimizer has moved outside the window
* **/
///three cases: 1.
if (info.x <= min.x) { // a new minimum; then write the old min
if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
min = info, min_pos = buf_pos;
} else if (buf_pos == min_pos) { // old min has moved outside the window
if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
///buf_pos == min_pos, means current minimizer has moved outside the window
///so for now we need to find a new minimizer at the current window (w k-mers)
for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers
if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
for (j = 0; j <= buf_pos; ++j)
if (min.x >= buf[j].x) min = buf[j], min_pos = j;
if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
for (j = 0; j <= buf_pos; ++j)
if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
}
}
if (++buf_pos == w) buf_pos = 0;
}
if (min.x != UINT64_MAX)
kv_push(ha_mz1_t, *p, min);
}