diff --git a/Assembly.cpp b/Assembly.cpp index 75d706f..976d304 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -848,7 +848,8 @@ void ha_overlap_and_correct(int round) b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1)); 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 + debug_adapter(&asm_opt, &R_INF); + 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); diff --git a/htab.cpp b/htab.cpp index dcbf9c9..fb40fe6 100644 --- a/htab.cpp +++ b/htab.cpp @@ -720,6 +720,48 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip return 0; } +void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs) +{ + int ret; + uint32_t i, m, pass, unpass; + gzFile fp = 0; + kseq_t *ks = NULL; + UC_Read ucr; + init_UC_Read(&ucr); + + for (i = m = pass = unpass = 0; i < (uint32_t)asm_opt->num_reads; ++i) + { + if ((fp = gzopen(asm_opt->read_file_names[i], "r")) == 0) continue; + ks = kseq_init(fp); + while ((ret = kseq_read(ks)) >= 0) + { + int l = ks->seq.l; + if((l - asm_opt->adapterLen*2) <= 0) continue; + recover_UC_Read(&ucr, rs, m); + fprintf(stderr, "l: %d, ucr.length: %lld, asm_opt->adapterLen: %d\n", + l, ucr.length, asm_opt->adapterLen); + if(memcmp(ucr.seq, ks->seq.s+asm_opt->adapterLen, ucr.length) == 0) + { + pass++; + } + else + { + unpass++; + } + m++; + } + kseq_destroy(ks); + gzclose(fp); + ks = NULL; + fp = 0; + } + + destory_UC_Read(&ucr); + + fprintf(stderr, "[M::%s::# reads: %u, # pass: %u, # unpass: %u\n]", __func__, m, pass, unpass); + exit(1); +} + static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt_t *p0, ha_ct_t *c0, const void *flt_tab, All_reads *rs, int64_t *n_seq) { ///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN|HAF_CREATE_NEW diff --git a/htab.h b/htab.h index c6feba6..613fbea 100644 --- a/htab.h +++ b/htab.h @@ -62,6 +62,7 @@ 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, kvec_t_u64_warp* dbg_ct); int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het); +void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs); static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function {