-z version number

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);

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.13-r308"
+#define HA_VERSION "0.13-r308-z"
 
 #define VERBOSE 0
 

htab.cpp

@@ -49,6 +49,7 @@ typedef struct {
 	int32_t pre;
 	int32_t n_thread;
 	int64_t chunk_size;
+	int adaLen;
 } yak_copt_t;
 
 void yak_copt_init(yak_copt_t *o)
@@ -565,159 +566,203 @@ static void worker_for_mz(void *data, long i, int tid)
 
 static void *worker_count(void *data, int step, void *in) // callback for kt_pipeline()
 {
-	pl_data_t *p = (pl_data_t*)data;
+    pl_data_t *p = (pl_data_t*)data;
-	if (step == 0) { // step 1: read a block of sequences
+    if (step == 0) { // step 1: read a block of sequences
-		int ret;
+        int ret;
-		st_data_t *s;
+        st_data_t *s;
-		CALLOC(s, 1);
+        CALLOC(s, 1);
-		s->p = p;
+        s->p = p;
-		s->n_seq0 = p->n_seq;
+        s->n_seq0 = p->n_seq;
-		if (p->rs_in && (p->flag & HAF_RS_READ)) {
+        if (p->rs_in && (p->flag & HAF_RS_READ)) {
-			while (p->n_seq < p->rs_in->total_reads) {
+            while (p->n_seq < p->rs_in->total_reads) {
-				if((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU)
+                if((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU)
-				{
+                {
-					++p->n_seq;
+                    ++p->n_seq;
-					continue;
+                    continue;
-				}
+                }
-				int l;
+                int l;
-				recover_UC_Read(&p->ucr, p->rs_in, p->n_seq);
+                recover_UC_Read(&p->ucr, p->rs_in, p->n_seq);
-				l = p->ucr.length;
+                l = p->ucr.length;
-				if (s->n_seq == s->m_seq) {
+                if (s->n_seq == s->m_seq) {
-					s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
+                    s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
-					REALLOC(s->len, s->m_seq);
+                    REALLOC(s->len, s->m_seq);
-					REALLOC(s->seq, s->m_seq);
+                    REALLOC(s->seq, s->m_seq);
-				}
+                }
-				MALLOC(s->seq[s->n_seq], l);
+                MALLOC(s->seq[s->n_seq], l);
-				memcpy(s->seq[s->n_seq], p->ucr.seq, l);
+                memcpy(s->seq[s->n_seq], p->ucr.seq, l);
-				s->len[s->n_seq++] = l;
+                s->len[s->n_seq++] = l;
-				++p->n_seq;
+                ++p->n_seq;
-				s->sum_len += l;
+                s->sum_len += l;
-				s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
+                s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
-				if (s->sum_len >= p->opt->chunk_size)
+                if (s->sum_len >= p->opt->chunk_size)
-					break;
+                    break;
-			}
+            }
-		} else {
+        } else {
-			while ((ret = kseq_read(p->ks)) >= 0) {
+            while ((ret = kseq_read(p->ks)) >= 0) {
-				int l = p->ks->seq.l;
+                int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen);
-				if (p->n_seq >= 1<<28) {
+                if(l <= 0) continue;
-					fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);
+
-					exit(1);
+                if (p->n_seq >= 1<<28) {
-				}
+                    fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);
-				if (p->rs_out) {
+                    exit(1);
-					///for 0-th count, just insert read length to R_INF, instead of read
+                }
-					if (p->flag & HAF_RS_WRITE_LEN) {
+                if (p->rs_out) {
-						assert(p->n_seq == p->rs_out->total_reads);
+                    ///for 0-th count, just insert read length to R_INF, instead of read
-						ha_insert_read_len(p->rs_out, l, p->ks->name.l);
+                    if (p->flag & HAF_RS_WRITE_LEN) {
-					} else if (p->flag & HAF_RS_WRITE_SEQ) {
+                        assert(p->n_seq == p->rs_out->total_reads);
-						int i, n_N;
+                        ha_insert_read_len(p->rs_out, l, p->ks->name.l);
-						assert(l == (int)p->rs_out->read_length[p->n_seq]);
+                    } else if (p->flag & HAF_RS_WRITE_SEQ) {
-						for (i = n_N = 0; i < l; ++i) // count number of ambiguous bases
+                        int i, n_N;
-							if (seq_nt4_table[(uint8_t)p->ks->seq.s[i]] >= 4)
+                        assert(l == (int)p->rs_out->read_length[p->n_seq]);
-								++n_N;
+                        for (i = n_N = 0; i < l; ++i) // count number of ambiguous bases
-						ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s, l, &p->rs_out->N_site[p->n_seq], n_N);
+                            if (seq_nt4_table[(uint8_t)p->ks->seq.s[i+p->opt->adaLen]] >= 4)
-						memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);
+                                ++n_N;
-					}
+                        ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N);
-				}
+                        memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);
-				///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);
+                ///for 0-th count, insert both seq and length to local block
-					REALLOC(s->len, s->m_seq);
+                if (s->n_seq == s->m_seq) {
-					REALLOC(s->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);
-				MALLOC(s->seq[s->n_seq], l);
+                    REALLOC(s->seq, s->m_seq);
-				memcpy(s->seq[s->n_seq], p->ks->seq.s, l);
+                }
-				s->len[s->n_seq++] = l;
+                MALLOC(s->seq[s->n_seq], l);
-				++p->n_seq;
+                memcpy(s->seq[s->n_seq], p->ks->seq.s+p->opt->adaLen, l);
-				s->sum_len += l;
+                s->len[s->n_seq++] = l;
-				s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
+                ++p->n_seq;
-				///p->opt->chunk_size is the block max size
+                s->sum_len += l;
-				if (s->sum_len >= p->opt->chunk_size)
+                s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
-					break;
+                ///p->opt->chunk_size is the block max size
-			}
+                if (s->sum_len >= p->opt->chunk_size)
-		}
+                    break;
-		if (s->sum_len == 0) free(s);
+            }
-		else return s;
+        }
-	} else if (step == 1) { // step 2: extract k-mers
+        if (s->sum_len == 0) free(s);
-		///s is the block of reads
+        else return s;
-		st_data_t *s = (st_data_t*)in;
+    } else if (step == 1) { // step 2: extract k-mers
-		///for 0-th counting, n_pre = 4096
+        ///s is the block of reads
-		int i, n_pre = 1<<p->opt->pre, m;
+        st_data_t *s = (st_data_t*)in;
-		// allocate the k-mer buffer
+        ///for 0-th counting, n_pre = 4096
-		CALLOC(s->buf, n_pre);
+        int i, n_pre = 1<<p->opt->pre, m;
-		m = (int)(s->nk * 1.2 / n_pre) + 1;
+        // allocate the k-mer buffer
-		//pre-allocate memory for each of 4096 buffer
+        CALLOC(s->buf, n_pre);
-		for (i = 0; i < n_pre; ++i) {
+        m = (int)(s->nk * 1.2 / n_pre) + 1;
-			s->buf[i].m = m;
+        //pre-allocate memory for each of 4096 buffer
-			///for 0-th counting, p->pt = NULL
+        for (i = 0; i < n_pre; ++i) {
-			if (p->pt) MALLOC(s->buf[i].b, m);
+            s->buf[i].m = m;
-			else MALLOC(s->buf[i].a, m);
+            ///for 0-th counting, p->pt = NULL
-		}
+            if (p->pt) MALLOC(s->buf[i].b, m);
-		// fill the buffer
+            else MALLOC(s->buf[i].a, m);
-		///for 0-th counting, p->opt->w == 1
+        }
-		if (p->opt->w == 1) { // enumerate all k-mers
+        // fill the buffer
-			///scan all reads
+        ///for 0-th counting, p->opt->w == 1
-			for (i = 0; i < s->n_seq; ++i) {
+        if (p->opt->w == 1) { // enumerate all k-mers
-				if (p->opt->is_HPC)
+            ///scan all reads
-					count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
+            for (i = 0; i < s->n_seq; ++i) {
-				else
+                if (p->opt->is_HPC)
-					count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
+                    count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
-				if (!p->is_store) free(s->seq[i]);
+                else
-			}
+                    count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
-		} else { // minimizers only
+                if (!p->is_store) free(s->seq[i]);
-			uint32_t j;
+            }
-			// compute minimizers
+        } else { // minimizers only
-			// s->n_seq is how many reads at this buffer
+            uint32_t j;
-			// s->mz && s->mz_buf are lists of minimzer vectors
+            // compute minimizers
-			CALLOC(s->mz, s->n_seq);
+            // s->n_seq is how many reads at this buffer
-			CALLOC(s->mz_buf, p->opt->n_thread);
+            // s->mz && s->mz_buf are lists of minimzer vectors
-			///calculate minimzers for each read, each read corresponds to one thread
+            CALLOC(s->mz, s->n_seq);
-			kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq);
+            CALLOC(s->mz_buf, p->opt->n_thread);
-			for (i = 0; i < p->opt->n_thread; ++i)
+            ///calculate minimzers for each read, each read corresponds to one thread
-				free(s->mz_buf[i].a);
+            kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq);
-			free(s->mz_buf);
+            for (i = 0; i < p->opt->n_thread; ++i)
-			// insert minimizers
+                free(s->mz_buf[i].a);
-			if (p->pt) {///insert whole minimizer
+            free(s->mz_buf);
-				for (i = 0; i < s->n_seq; ++i)
+            // insert minimizers
-					for (j = 0; j < s->mz[i].n; ++j)
+            if (p->pt) {///insert whole minimizer
-						pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]);
+                for (i = 0; i < s->n_seq; ++i)
-			} else {///just insert the hash key of minimizer
+                    for (j = 0; j < s->mz[i].n; ++j)
-				for (i = 0; i < s->n_seq; ++i)
+                        pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]);
-					for (j = 0; j < s->mz[i].n; ++j)
+            } else {///just insert the hash key of minimizer
-						ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
+                for (i = 0; i < s->n_seq; ++i)
-			}
+                    for (j = 0; j < s->mz[i].n; ++j)
-			for (i = 0; i < s->n_seq; ++i) {
+                        ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
-				free(s->mz[i].a);
+            }
-				if (!p->is_store) free(s->seq[i]);
+            for (i = 0; i < s->n_seq; ++i) {
-			}
+                free(s->mz[i].a);
-			free(s->mz);
+                if (!p->is_store) free(s->seq[i]);
-		}
+            }
-		///just clean seq
+            free(s->mz);
-		free(s->seq); free(s->len);
+        }
-		s->seq = 0, s->len = 0;
+        ///just clean seq
-		return s;
+        free(s->seq); free(s->len);
-	} else if (step == 2) { // step 3: insert k-mers to hash table
+        s->seq = 0, s->len = 0;
-		st_data_t *s = (st_data_t*)in;
+        return s;
-		int i, n = 1<<p->opt->pre;
+    } else if (step == 2) { // step 3: insert k-mers to hash table
-		uint64_t n_ins = 0;
+        st_data_t *s = (st_data_t*)in;
-		///for 0-th counting, p->pt = NULL
+        int i, n = 1<<p->opt->pre;
-		kt_for(p->opt->n_thread, worker_for_insert, s, n);
+        uint64_t n_ins = 0;
-		///n_ins is number of distinct k-mers
+        ///for 0-th counting, p->pt = NULL
-		for (i = 0; i < n; ++i) {
+        kt_for(p->opt->n_thread, worker_for_insert, s, n);
-			n_ins += s->buf[i].n_ins;
+        ///n_ins is number of distinct k-mers
-			if (p->pt) free(s->buf[i].b);
+        for (i = 0; i < n; ++i) {
-			else free(s->buf[i].a);
+            n_ins += s->buf[i].n_ins;
-		}
+            if (p->pt) free(s->buf[i].b);
-		if (p->ct) p->ct->tot += n_ins;
+            else free(s->buf[i].a);
-		if (p->pt) p->pt->tot_pos += n_ins;
+        }
-		free(s->buf);
+        if (p->ct) p->ct->tot += n_ins;
-		#if 0
+        if (p->pt) p->pt->tot_pos += n_ins;
-		fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__,
+        free(s->buf);
-				yak_realtime(), yak_cpu_usage(), (long)s->n_seq0 + s->n_seq,
+        #if 0
-				(long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers");
+        fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__,
-		#endif
+                yak_realtime(), yak_cpu_usage(), (long)s->n_seq0 + s->n_seq,
-		free(s);
+                (long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers");
-	}
+        #endif
-	return 0;
+        free(s);
+    }
+    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)
@@ -769,34 +814,35 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 
 ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const void *flt_tab, All_reads *rs)
 {
-	int i;
+    int i;
-	int64_t n_seq = 0;
+    int64_t n_seq = 0;
-	yak_copt_t opt;
+    yak_copt_t opt;
-	ha_ct_t *h = 0;
+    ha_ct_t *h = 0;
-	assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both
+    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
+    ///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN
-	if (rs) {
+    if (rs) {
-		if (flag & HAF_RS_WRITE_LEN)
+        if (flag & HAF_RS_WRITE_LEN)
-			init_All_reads(rs);
+            init_All_reads(rs);
-		else if (flag & HAF_RS_WRITE_SEQ)
+        else if (flag & HAF_RS_WRITE_SEQ)
-			malloc_All_reads(rs);
+            malloc_All_reads(rs);
-	}
+    }
-	yak_copt_init(&opt);
+    yak_copt_init(&opt);
-	opt.k = asm_opt->k_mer_length;
+    opt.k = asm_opt->k_mer_length;
-	///always 0
+    ///always 0
-	opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC);
+    opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC);
-	///for ft-counting, shoud be 1
+    ///for ft-counting, shoud be 1
-	opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win;
+    opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win;
-	///for ft-counting, shoud be 37
+    ///for ft-counting, shoud be 37
-	///for ha_pt_gen, shoud be 0
+    ///for ha_pt_gen, shoud be 0
-	opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift;
+    opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift;
-	opt.n_thread = asm_opt->thread_num;
+    opt.n_thread = asm_opt->thread_num;
-	///asm_opt->num_reads is the number of fastq files
+    opt.adaLen = asm_opt->adapterLen;
-	for (i = 0; i < asm_opt->num_reads; ++i)
+    ///asm_opt->num_reads is the number of fastq files
-		h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, &n_seq);
+    for (i = 0; i < asm_opt->num_reads; ++i)
-	if (h && opt.bf_shift > 0)
+        h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, &n_seq);
-		ha_ct_destroy_bf(h);
+    if (h && opt.bf_shift > 0)
-	return h;
+        ha_ct_destroy_bf(h);
+    return h;
 }
 
 /***************************

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
 {