r254: added --ex-list and --ex-iter

Assembly.cpp

@@ -1159,10 +1159,15 @@ void ha_overlap_final(void)
 
 int ha_assemble(void)
 {
+	extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
 	int r, hom_cov = -1, ovlp_loaded = 0;
 	if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
 		ovlp_loaded = 1;
 		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
+		if (asm_opt.extract_list) {
+			ha_extract_print_list(&R_INF, asm_opt.extract_iter, asm_opt.extract_list);
+			exit(0);
+		}
 		if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt);
         ///if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt);
 		if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();

CommandLines.cpp

@@ -17,6 +17,8 @@ static ko_longopt_t long_options[] = {
 	{ "skip-triobin",  ko_no_argument, 304 },
 	{ "max-od-ec",     ko_no_argument, 305 },
 	{ "max-od-final",  ko_no_argument, 306 },
+	{ "ex-list",       ko_required_argument, 307 },
+	{ "ex-iter",       ko_required_argument, 308 },
 	{ 0, 0, 0 }
 };
 
@@ -400,6 +402,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 		else if (c == 304) asm_opt->flag |= HA_F_SKIP_TRIOBIN;
 		else if (c == 305) asm_opt->max_ov_diff_ec = atof(opt.arg);
 		else if (c == 306) asm_opt->max_ov_diff_final = atof(opt.arg);
+		else if (c == 307) asm_opt->extract_list = opt.arg;
+		else if (c == 308) asm_opt->extract_iter = atoi(opt.arg);
         else if (c == 'l')
         {   ///0: disable purge_dup; 1: purge containment; 2: purge overlap
             asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.7-dirty-r253"
+#define HA_VERSION "0.7-dirty-r254"
 
 #define VERBOSE 0
 
@@ -26,6 +26,8 @@ typedef struct {
     char* required_read_name;
 	char *fn_bin_yak[2];
 	char *fn_bin_list[2];
+	char *extract_list;
+	int extract_iter;
     int thread_num;
     int k_mer_length;
 	int mz_win;

Makefile

@@ -4,7 +4,7 @@ CPPFLAGS=
 INCLUDES=
 OBJS=		CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
 			POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \
-			htab.o hist.o sketch.o anchor.o sys.o
+			htab.o hist.o sketch.o anchor.o extract.o sys.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm
 
@@ -51,10 +51,12 @@ Process_Read.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
 Purge_Dups.o: ksort.h Purge_Dups.h kvec.h kdq.h Overlaps.h Hash_Table.h
 Purge_Dups.o: htab.h Process_Read.h CommandLines.h Correct.h
 Purge_Dups.o: Levenshtein_distance.h POA.h kthread.h
-Trio.o: khashl.h kthread.h Process_Read.h Overlaps.h kvec.h kdq.h
+Trio.o: khashl.h kthread.h kseq.h Process_Read.h Overlaps.h kvec.h kdq.h
 Trio.o: CommandLines.h htab.h
 anchor.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
 anchor.o: ksort.h Hash_Table.h
+extract.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h khashl.h
+extract.o: kseq.h
 hist.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
 htab.o: kthread.h khashl.h kseq.h ksort.h htab.h Process_Read.h Overlaps.h
 htab.o: kvec.h kdq.h CommandLines.h

Overlaps.cpp

@@ -10506,7 +10506,7 @@ long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negativ
                 {
                     to_replace = 1;
                 }
-                else if((cur_weight == max_weight))
+                else if (cur_weight == max_weight)
                 {
                     if(ll > base_maxLen)
                     {

extract.cpp

@@ -0,0 +1,175 @@
+#include <zlib.h>
+#include <string.h>
+#include "Process_Read.h"
+#include "khashl.h"
+#include "kseq.h"
+
+typedef const char *cstr_t;
+KHASHL_CMAP_INIT(KH_LOCAL, strmap_t, ss, cstr_t, int, kh_hash_str, kh_eq_str)
+KSTREAM_INIT(gzFile, gzread, 65536)
+
+#define GFA_MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr))))
+#define GFA_REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr))))
+
+char *gfa_strdup(const char *src)
+{
+	int32_t len;
+	char *dst;
+	len = strlen(src);
+	GFA_MALLOC(dst, len + 1);
+	memcpy(dst, src, len + 1);
+	return dst;
+}
+
+char *gfa_strndup(const char *src, size_t n)
+{
+	char *dst;
+	GFA_MALLOC(dst, n + 1);
+	strncpy(dst, src, n);
+	dst[n] = 0;
+	return dst;
+}
+
+char **gv_read_list(const char *o, int *n_)
+{
+	int n = 0, m = 0;
+	char **s = 0;
+	*n_ = 0;
+	if (*o != '@') {
+		const char *q = o, *p;
+		for (p = q;; ++p) {
+			if (*p == ',' || *p == 0) {
+				if (n == m) {
+					m = m? m<<1 : 16;
+					GFA_REALLOC(s, m);
+				}
+				s[n++] = gfa_strndup(q, p - q);
+				if (*p == 0) break;
+				q = p + 1;
+			}
+		}
+	} else {
+		gzFile fp;
+		kstream_t *ks;
+		kstring_t str = {0,0,0};
+		int dret;
+
+		fp = gzopen(o + 1, "r");
+		if (fp == 0) return 0;
+		ks = ks_init(fp);
+		while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) {
+			char *p;
+			for (p = str.s; *p && !isspace(*p); ++p);
+			if (n == m) {
+				m = m? m<<1 : 16;
+				GFA_REALLOC(s, m);
+			}
+			s[n++] = gfa_strndup(str.s, p - str.s);
+		}
+		ks_destroy(ks);
+		gzclose(fp);
+	}
+	if (s) s = (char**)realloc(s, n * sizeof(char*));
+	*n_ = n;
+	return s;
+}
+
+void ha_extract_print(const All_reads *rs, int n_rounds, int n, char **list)
+{
+	strmap_t *h;
+	khint_t k;
+	int i, absent, m, l, max_len = 0;
+	uint64_t j;
+	char *s = 0;
+	const ma_hit_t_alloc *ov[2] = { rs->paf, rs->reverse_paf };
+	FILE *fp = stdout;
+
+	for (j = 0; j < rs->total_reads; ++j) {
+		if (max_len < (int)Get_NAME_LENGTH(*rs, j))
+			max_len = Get_NAME_LENGTH(*rs, j);
+	}
+	GFA_MALLOC(s, max_len + 1);
+
+	h = ss_init();
+	for (i = 0; i < n; ++i) {
+		k = ss_put(h, gfa_strdup(list[i]), &absent);
+		kh_val(h, k) = -1;
+	}
+
+	for (m = 0; m < n_rounds; ++m) {
+		for (j = 0; j < rs->total_reads; ++j) {
+			for (l = 0; l < 2; ++l) {
+				const ma_hit_t_alloc *o = &ov[l][j];
+				for (i = 0; i < (int)o->length; ++i) {
+					uint64_t q = Get_qn(o->buffer[i]);
+					uint64_t t = Get_tn(o->buffer[i]);
+					int q_hit = 0, t_hit = 0;
+					strncpy(s, Get_NAME(*rs, q), Get_NAME_LENGTH(*rs, q)); s[Get_NAME_LENGTH(*rs, q)] = 0;
+					k = ss_get(h, s);
+					q_hit = (k < kh_end(h) && kh_val(h, k) < m);
+					strncpy(s, Get_NAME(*rs, t), Get_NAME_LENGTH(*rs, t)); s[Get_NAME_LENGTH(*rs, t)] = 0;
+					k = ss_get(h, s);
+					t_hit = (k < kh_end(h) && kh_val(h, k) < m);
+					if ((!q_hit && !t_hit) || (q_hit && t_hit)) continue;
+					if (!q_hit) {
+						char *tmp = gfa_strndup(Get_NAME(*rs, q), Get_NAME_LENGTH(*rs, q));
+						k = ss_put(h, tmp, &absent);
+						if (absent) kh_val(h, k) = m;
+						else free(tmp);
+					}
+					if (!t_hit) {
+						char *tmp = gfa_strndup(Get_NAME(*rs, t), Get_NAME_LENGTH(*rs, t));
+						k = ss_put(h, tmp, &absent);
+						if (absent) kh_val(h, k) = m;
+						else free(tmp);
+					}
+				}
+			}
+		}
+	}
+
+	for (j = 0; j < rs->total_reads; ++j) {
+		for (l = 0; l < 2; ++l) {
+			const ma_hit_t_alloc *o = &ov[l][j];
+			for (i = 0; i < (int)o->length; ++i) {
+				uint64_t q = Get_qn(o->buffer[i]);
+				uint64_t t = Get_tn(o->buffer[i]);
+				int q_hit = 0, t_hit = 0;
+				strncpy(s, Get_NAME(*rs, q), Get_NAME_LENGTH(*rs, q)); s[Get_NAME_LENGTH(*rs, q)] = 0;
+				k = ss_get(h, s);
+				q_hit = (k < kh_end(h));
+				strncpy(s, Get_NAME(*rs, t), Get_NAME_LENGTH(*rs, t)); s[Get_NAME_LENGTH(*rs, t)] = 0;
+				k = ss_get(h, s);
+				t_hit = (k < kh_end(h));
+				if (!q_hit && !t_hit) continue;
+				fwrite(Get_NAME(*rs, q), 1, Get_NAME_LENGTH(*rs, q), fp);
+				fwrite("\t", 1, 1, fp);
+				fprintf(fp, "%lu\t", (unsigned long)Get_READ_LENGTH(*rs, q));
+				fprintf(fp, "%d\t", Get_qs(o->buffer[i]));
+				fprintf(fp, "%d\t", Get_qe(o->buffer[i]));
+				fputs(o->buffer[i].rev? "-\t" : "+\t", fp);
+				fwrite(Get_NAME(*rs, t), 1, Get_NAME_LENGTH(*rs, t), fp);
+				fwrite("\t", 1, 1, fp);
+				fprintf(fp, "%lu\t", (unsigned long)Get_READ_LENGTH(*rs, t));
+				fprintf(fp, "%d\t", Get_ts(o->buffer[i]));
+				fprintf(fp, "%d\t%d\t%d\t%d\n", Get_te(o->buffer[i]), o->buffer[i].ml, o->buffer[i].bl, !l);
+			}
+		}
+	}
+
+	for (k = 0; k != kh_end(h); ++k)
+		if (kh_exist(h, k))
+			free((char*)kh_key(h, k));
+	ss_destroy(h);
+	free(s);
+}
+
+void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o)
+{
+	int i, n;
+	char **list;
+	list = gv_read_list(o, &n);
+	ha_extract_print(rs, n_rounds, n, list);
+	for (i = 0; i < n; ++i) free(list[i]);
+	free(list);
+}