diff --git a/Assembly.cpp b/Assembly.cpp index b56c29c..d2b9134 100644 --- a/Assembly.cpp +++ b/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(); diff --git a/CommandLines.cpp b/CommandLines.cpp index f39bc98..b165697 100644 --- a/CommandLines.cpp +++ b/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); diff --git a/CommandLines.h b/CommandLines.h index 72aca9d..b453e6d 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#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; diff --git a/Makefile b/Makefile index 1cb14d0..250c1ae 100644 --- a/Makefile +++ b/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 diff --git a/Overlaps.cpp b/Overlaps.cpp index b27c0de..6ff6259 100644 --- a/Overlaps.cpp +++ b/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) { diff --git a/extract.cpp b/extract.cpp new file mode 100644 index 0000000..d139910 --- /dev/null +++ b/extract.cpp @@ -0,0 +1,175 @@ +#include +#include +#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); +}