update r312

Added citation

CommandLines.cpp

@@ -30,6 +30,8 @@ static ko_longopt_t long_options[] = {
     { "h2",      ko_required_argument, 315 },
     { "enzyme",      ko_required_argument, 316 },
     { "b-cov",      ko_required_argument, 317 },
+    { "h-cov",      ko_required_argument, 318 },
+    { "m-rate",      ko_required_argument, 319 },
 	{ 0, 0, 0 }
 };
 
@@ -68,7 +70,12 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "    --lowQ      INT\n");
     fprintf(stderr, "                output contig regions with >=INT%% inconsistency in BED format; 0 to disable [%d]\n", asm_opt->bed_inconsist_rate);
     fprintf(stderr, "    --b-cov     INT\n");
-    fprintf(stderr, "                break contigs at breakpoints with coverage drop at <INT-fold coverage [%d]\n", asm_opt->break_cov);
+    fprintf(stderr, "                break contigs at positions with <INT-fold coverage; work with '--m-rate'; 0 to disable [%d]\n", asm_opt->b_low_cov);
+    fprintf(stderr, "    --h-cov     INT\n");
+    fprintf(stderr, "                break contigs at positions with >INT-fold coverage; work with '--m-rate'; -1 to disable [%d]\n", asm_opt->b_high_cov);
+    fprintf(stderr, "    --m-rate    FLOAT\n");
+    fprintf(stderr, "                break contigs at positions with <=FLOAT*coverage exact overlaps;\n");
+    fprintf(stderr, "                only work with '--b-cov' or '--h-cov'[%.2f]\n", asm_opt->m_rate);
     
 //	fprintf(stderr, "    --pri-range INT1[,INT2]\n");
 //	fprintf(stderr, "                keep contigs with coverage in this range in p_ctg.gfa; -1 to disable [auto,inf]\n");
@@ -153,7 +160,9 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->hic_inconsist_rate = 30;
     ///asm_opt->bub_mer_length = 3;
     asm_opt->bub_mer_length = 1000000;
-    asm_opt->break_cov = 0;
+    asm_opt->b_low_cov = 0;
+    asm_opt->b_high_cov = -1;
+    asm_opt->m_rate = 0.75;
 }
 
 void destory_enzyme(enzyme* f)
@@ -419,6 +428,30 @@ int check_option(hifiasm_opt_t* asm_opt)
         return 0;
     }
 
+    if(asm_opt->b_low_cov < 0)
+    {
+        fprintf(stderr, "[ERROR] must >= 0 (--b-cov)\n");
+        return 0;
+    }
+
+    if(asm_opt->b_high_cov != -1 && asm_opt->b_high_cov < 0)
+    {
+        fprintf(stderr, "[ERROR] must >= 0 (--h-cov)\n");
+        return 0;
+    }
+
+    if(asm_opt->m_rate < 0)
+    {
+        fprintf(stderr, "[ERROR] must >= 0 (--m-rate)\n");
+        return 0;
+    }
+
+    if(asm_opt->b_high_cov != -1 && asm_opt->b_high_cov <= asm_opt->b_low_cov)
+    {
+        fprintf(stderr, "[ERROR] [--h-cov] must >= [--b-cov]\n");
+        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);
@@ -592,7 +625,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 314) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[0]), 0);
         else if (c == 315) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[1]), 0);
         else if (c == 316) get_hic_enzymes(opt.arg, &(asm_opt->hic_enzymes), 1);
-        else if (c == 317) asm_opt->break_cov = atoi(opt.arg);
+        else if (c == 317) asm_opt->b_low_cov = atoi(opt.arg);
+        else if (c == 318) asm_opt->b_high_cov = atoi(opt.arg);
+        else if (c == 319) asm_opt->m_rate = atof(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.14-r309"
+#define HA_VERSION "0.14-r312"
 
 #define VERBOSE 0
 
@@ -49,7 +49,9 @@ typedef struct {
 	double max_ov_diff_final;
 	int hom_cov;
     int het_cov;
-    int break_cov;
+    int b_low_cov;
+    int b_high_cov;
+    double m_rate;
 	int max_n_chain; // fall-back max number of chains to consider
 	int min_hist_kmer_cnt;
     int load_index_from_disk;

Makefile

@@ -43,7 +43,7 @@ Assembly.o: kthread.h
 CommandLines.o: CommandLines.h ketopt.h
 Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h
 Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h
-Correct.o: ksw2.h
+Correct.o: ksw2.h ksort.h
 Hash_Table.o: Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h
 Hash_Table.o: CommandLines.h ksort.h
 Levenshtein_distance.o: Levenshtein_distance.h

Overlaps.h

@@ -477,7 +477,7 @@ void destory_R_to_U(R_to_U* x);
 void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8_t* flag);
 void get_R_to_U(R_to_U* x, uint32_t rID, uint32_t* uID, uint32_t* is_Unitig);
 void transfor_R_to_U(R_to_U* x);
-void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, int require_equal_nv, int test_tangle);
+void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, kvec_asg_arc_t_warp* edge, int require_equal_nv, int test_tangle);
 int asg_pop_bubble_primary(asg_t *g, int max_dist);
 long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *g, long long circleLen, int is_drop);
 

Purge_Dups.cpp

@@ -1493,11 +1493,10 @@ void quick_LIS(asg_arc_t_offset* x, uint32_t n, kvec_t_i32_warp* tailIndex, kvec
     tailIndex->a.n = prevIndex->a.n = 0;
     if(n == 0) return;
 
-    
     kv_resize(int32_t, tailIndex->a, n);
     kv_resize(int32_t, prevIndex->a, n);
 
-    long long len = 1, i, pos, m;
+    long long len = 1, i, pos, m; ///the length of chain must be >=1
     tailIndex->a.a[0] = 0;
     prevIndex->a.a[0] = -1;
 
@@ -1521,8 +1520,7 @@ void quick_LIS(asg_arc_t_offset* x, uint32_t n, kvec_t_i32_warp* tailIndex, kvec
             // future subsequence 
             // It will replace ceil value in tailIndices 
             pos = GetCeilIndex(x, tailIndex, -1, len - 1, Get_yOff(x[i].Off)); 
-  
-            prevIndex->a.a[i] = tailIndex->a.a[pos - 1]; 
+            prevIndex->a.a[i] = pos > 0? tailIndex->a.a[pos - 1] : -1; 
             tailIndex->a.a[pos] = i; 
         } 
     }
@@ -1627,7 +1625,6 @@ uint32_t* xBeg, uint32_t* xEnd, uint32_t* yBeg, uint32_t* yEnd)
     u_buffer->a.n = m;
 
     ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "after");
-
     quick_LIS(u_buffer->a.a, u_buffer->a.n, tailIndex, prevIndex);
 
     if(tailIndex->a.n == 0) return;
@@ -1931,6 +1928,7 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon
         ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, min_ovlp, 
         xUid, yUid, u_buffer, tailIndex, prevIndex, r_x_pos_beg, r_x_pos_end, r_y_pos_beg, 
         r_y_pos_end);
+
         if(hap_can->index_end == XCY && yReads->len > (xReads->len*2)) return NON_PLOID;
         if(hap_can->index_end == YCX && xReads->len > (yReads->len*2)) return NON_PLOID;
         if(hap_can->index_end == (uint32_t)-1) return NON_PLOID;
@@ -2979,6 +2977,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid)
             {
                 continue;
             }
+
             
             if(Get_match(hap_can) < Get_match(u_can->a.a[k]))
             {
@@ -4136,6 +4135,7 @@ long long get_specific_contig_length(asg_t *g, uint8_t *del)
         totalLen += len;
     }
     //kdq_destroy(uint64_t, q);
+    free(mark);
     return totalLen;
 }
 
@@ -4181,6 +4181,7 @@ uint32_t minLen, double purge_threshold)
     asg_t* nsg = ug->g;
     uint64_t v, k, total_bases = 0, alter_bases = 0, primary_bases = 0, purge_bases = 0;
     kt_for(asm_opt.thread_num, hap_alignment_advance_worker, hap_buf, nsg->n_seq);
+    
     filter_hap_overlaps_by_length(all_ovlp, minLen);
     normalize_hap_overlaps_advance(all_ovlp, back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
 
@@ -4295,7 +4296,7 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link)
     init_hap_alignment_struct_pip(&hap_buf, asm_opt.thread_num, nsg->n_seq, ug, read_g,  
     sources, reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp, 
     0.05, &all_ovlp);
-
+    
     if(hap_buf.cov_threshold < 0)
     {
         if(if_ploid_sample(ug, read_g, ruIndex, sources, reverse_sources, coverage_cut, 

README.md

@@ -192,4 +192,8 @@ page](https://github.com/chhylp123/hifiasm/issues).
 
 ## Limitations
 
-1. Purging haplotig duplications may introduce misassemblies.
+1. Purging haplotig duplications may introduce misassemblies.
+
+## Citation
+
+Cheng, H., Concepcion, G.T., Feng, X., Zhang, H., Li H. Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. Nat Methods 18, 170–175 (2021). https://doi.org/10.1038/s41592-020-01056-5

hic.cpp

@@ -444,7 +444,7 @@ void hc_pt_t_gen_single(hc_pt1_t* pt, uint64_t* up_bound)
 int write_hc_pt_index(ha_ug_index* idx, char* file_name)
 {
     char* gfa_name = (char*)malloc(strlen(file_name)+25);
-    sprintf(gfa_name, "%s.hc_tlb", file_name);
+    sprintf(gfa_name, "%s.hic.tlb.bin", file_name);
     FILE* fp = fopen(gfa_name, "w");
     if (!fp) {
         free(gfa_name);
@@ -481,7 +481,7 @@ int load_hc_pt_index(ha_ug_index** r_idx, char* file_name)
     uint64_t flag = 0;
     double index_time = yak_realtime();
     char* gfa_name = (char*)malloc(strlen(file_name)+25);
-    sprintf(gfa_name, "%s.hc_tlb", file_name);
+    sprintf(gfa_name, "%s.hic.tlb.bin", file_name);
     FILE* fp = fopen(gfa_name, "r");
     if (!fp) {
         free(gfa_name);
@@ -2142,6 +2142,7 @@ void destory_pdq(pdq* q)
 {
     kv_destroy(q->x);
     kv_destroy(q->dis);
+    kv_destroy(q->vis);
 }
 
 void reset_pdq(pdq* q)
@@ -2889,7 +2890,7 @@ int load_hc_links(hc_links* link, const char *fn)
 void write_hc_hits(kvec_pe_hit* hits, const char *fn)
 {
     char *buf = (char*)calloc(strlen(fn) + 25, 1);
-	sprintf(buf, "%s.hic.lk", fn);
+	sprintf(buf, "%s.hic.lk.bin", fn);
     FILE* fp = fopen(buf, "w");
 
     fwrite(&hits->a.n, sizeof(hits->a.n), 1, fp);
@@ -2903,7 +2904,7 @@ int load_hc_hits(kvec_pe_hit* hits, const char *fn)
 {
     uint64_t flag = 0;
     char *buf = (char*)calloc(strlen(fn) + 25, 1);
-	sprintf(buf, "%s.hic.lk", fn);
+	sprintf(buf, "%s.hic.lk.bin", fn);
 
     FILE* fp = NULL; 
     fp = fopen(buf, "r"); 
@@ -7369,7 +7370,7 @@ void init_chain_hic_warp(ma_ug_t* ug, hc_links* link, bubble_type* bub, chain_hi
         }
     }
 
-    fprintf(stderr, "# chain: %u, # pre chain: %u\n", m, (uint32_t)(*c_w).n);
+    ///fprintf(stderr, "# chain: %u, # pre chain: %u\n", m, (uint32_t)(*c_w).n);
     (*c_w).n = m;
     for (i = 0; i < (*c_w).n; i++)
     {
@@ -11272,7 +11273,7 @@ int alignment_worker_pipeline(sldat_t* sl, const enzyme *fn1, const enzyme *fn2)
 
         kt_pipeline(3, worker_pipeline, sl, 3);
 
-        fprintf(stderr, "fn1->a[i]: %s, fn2->a[i]: %s, sl->hits.a.n: %u\n", fn1->a[i], fn2->a[i], (uint32_t)sl->hits.a.n);
+        ///fprintf(stderr, "fn1->a[i]: %s, fn2->a[i]: %s, sl->hits.a.n: %u\n", fn1->a[i], fn2->a[i], (uint32_t)sl->hits.a.n);
         
         
         kseq_destroy(sl->ks1);
@@ -11281,11 +11282,11 @@ int alignment_worker_pipeline(sldat_t* sl, const enzyme *fn1, const enzyme *fn2)
         gzclose(fp2);
     }
 
-    fprintf(stderr, "+sl->hits.a.n: %u\n", (uint32_t)sl->hits.a.n);
+    ///fprintf(stderr, "+sl->hits.a.n: %u\n", (uint32_t)sl->hits.a.n);
 
     dedup_hits(&(sl->hits));
 
-    fprintf(stderr, "-sl->hits.a.n: %u\n", (uint32_t)sl->hits.a.n);
+    ///fprintf(stderr, "-sl->hits.a.n: %u\n", (uint32_t)sl->hits.a.n);
     
     return 1;
 }
@@ -11401,127 +11402,6 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx)
     return 1;
 }
 
-int hic_short_align_back(const char *fn1, const char *fn2, ha_ug_index* idx)
-{
-    double index_time = yak_realtime();
-    sldat_t sl;
-    gzFile fp1, fp2;
-    kvec_hc_edge back_hc_edge;
-    kv_init(back_hc_edge.a);
-    if ((fp1 = gzopen(fn1, "r")) == 0) return 0;
-    if ((fp2 = gzopen(fn2, "r")) == 0) return 0;
-    sl.ks1 = kseq_init(fp1);
-    sl.ks2 = kseq_init(fp2);
-    sl.idx = idx;
-    sl.link = idx->link;
-    sl.chunk_size = 20000000;
-    sl.n_thread = asm_opt.thread_num;
-    sl.total_base = sl.total_pair = 0;
-    idx->max_cnt = 5;
-    kv_init(sl.hits.a);
-    
-    if(!load_hc_hits(&sl.hits, asm_opt.output_file_name))
-    {
-        /*******************************for debug************************************/
-        // load_reads(&R1, fn1);
-        // test_reads(&R1, fn1);
-        // load_reads(&R2, fn2);
-        // test_reads(&R1, fn1);
-        /*******************************for debug************************************/
-        fprintf(stderr, "+sl->hits.a.n: %u\n", (uint32_t)sl.hits.a.n);
-        kt_pipeline(3, worker_pipeline, &sl, 3);
-        fprintf(stderr, "-sl->hits.a.n: %u\n", (uint32_t)sl.hits.a.n);
-        fprintf(stderr, "fn1: %s, fn2: %s\n", fn1, fn2);
-
-        /*******************************for debug************************************/
-        // sort_hits(&sl.hits);
-        // print_hits(idx, &sl.hits, fn1);
-        /*******************************for debug************************************/
-        dedup_hits(&sl.hits);
-        write_hc_hits(&sl.hits, asm_opt.output_file_name);
-    }
-
-    fprintf(stderr, "u.n: %d, uID_bits: %lu, pos_bits: %lu, sl.hits.a.n: %u\n", (uint32_t)idx->ug->u.n, idx->uID_bits, idx->pos_bits, (uint32_t)sl.hits.a.n);
-
-    H_partition hap;
-    MT M;
-    init_MT(&M, idx->ug->g->n_seq<<1);
-    bubble_type bub; 
-    memset(&bub, 0, sizeof(bubble_type));
-    bub.round_id = 0; bub.n_round = 2;
-    for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++)
-    {
-        identify_bubbles(idx->ug, &bub, idx->link);
-        if(bub.round_id == 0)
-        {
-            collect_hc_links(sl.idx, &sl.hits, idx->link, &bub, &M);
-            collect_hc_reverse_links(idx->link, idx->ug, &bub);
-        }
-        init_hic_p((ha_ug_index*)sl.idx, &sl.hits, idx->link, &bub, &back_hc_edge, &M, &hap, 0);
-        ///init_hic_p_new((ha_ug_index*)sl.idx, &sl.hits, idx->link, &bub, &back_hc_edge, &M);
-        reset_H_partition(&hap, (bub.round_id == 0? 1 : 0));
-        init_contig_partition(&hap, idx, &bub);
-        phasing_improvement(&hap, &(hap.g_p), idx, &bub);
-        label_unitigs(&(hap.g_p), idx->ug);
-
-        ///print_hc_links(idx->link, 0, &hap);
-    }
-
-    cluster_contigs(&bub, idx, &sl.hits, &M, &hap);
-
-    destory_MT(&M);
-
-    ///print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, idx->link, idx);
-    ///print_hits(idx, &sl.hits, fn1);
-    
-
-    ///print_debug_bubble_graph(&bub, idx->ug, asm_opt.output_file_name);
-    // print_bubble_chain(&bub);
-    // print_hc_links(idx->link, 0, &hap);
-    
-    ///print_contig_partition(&hap, "final");
-
-    // uint32_t i;
-    // for (i = 0; i < idx->ug->g->n_seq; i++)
-    // {
-    //     fprintf(stderr, "utg%.6ul, index: %u\n", (int)(i+1), bub.index[i]);
-    // }
-    
-
-
-
-    destory_contig_partition(&hap);
-    kv_destroy(back_hc_edge.a);
-    return 1;
-    
-    /*******************************for debug************************************/
-    // destory_reads(&R1);
-    // destory_reads(&R2);
-    /*******************************for debug************************************/
-    print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, idx->link, idx);
-    collect_hc_reverse_links(idx->link, idx->ug, &bub);
-    normalize_hc_links(idx->link);
-    /*******************************for debug************************************/
-    ///print_hc_links(&link);
-    /*******************************for debug************************************/
-    min_cut_t* cut = clean_hap(idx->link, &bub, idx->ug);
-    ///print_bubbles(idx->ug, &bub, NULL, &link, idx);
-    G_partition* gp = clean_bubbles(idx->link, &bub, cut, idx->ug);
-    ///print_hc_links(&link);
-
-    destory_min_cut_t(cut); free(cut);
-    destory_G_partition(gp); free(gp);
-    kv_destroy(sl.hits.a);
-    destory_bubbles(&bub);
-    kseq_destroy(sl.ks1);
-    kseq_destroy(sl.ks2);
-    gzclose(fp1);
-    gzclose(fp2);
-    
-    fprintf(stderr, "[M::%s::%.3f] processed %lu pairs; %lu bases\n", __func__, yak_realtime()-index_time, sl.total_pair, sl.total_base);
-    return 1;
-}
-
 
 void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link)
 {
@@ -11533,8 +11413,7 @@ void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link)
     ug_index->read_g = read_g;
     ug_index->link = link;
     ///test_unitig_index(ug_index, ug);
-    ///hic_short_align(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index);
-    hic_short_align_back(asm_opt.hic_reads[0]->a[0], asm_opt.hic_reads[1]->a[0], ug_index);
+    hic_short_align(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index);
     
     destory_hc_pt_index(ug_index);
 }

hifiasm.1

@@ -1,4 +1,4 @@
-.TH hifiasm 1 "19 July 2020" "hifiasm-0.9 (r289)" "Bioinformatics tools"
+.TH hifiasm 1 "13 Feb 2021" "hifiasm-0.14 (r310)" "Bioinformatics tools"
 
 .SH NAME
 .PP
@@ -212,6 +212,32 @@ with suffix
 .B lowQ.bed
 [70]. Set 0 to disable. 
 
+
+.TP
+.BI --b-cov \ INT
+Break contigs at potential misassemblies with <INT-fold coverage [0].
+Work with 
+.B --m-rate.
+Set 0 to disable. 
+
+.TP
+.BI --h-cov \ INT
+Break contigs at potential misassemblies with >INT-fold coverage [-1].
+Work with 
+.B --m-rate.
+Set -1 to disable. 
+
+.TP
+.BI --m-rate \ FLOAT
+Break contigs with <=FLOAT*coverage exact overlaps [0.75].
+Only work with 
+.B --b-cov
+and
+.B --h-cov.
+
+
+
+
 .SS Trio-partition options
 
 .TP 10
@@ -289,6 +315,17 @@ For ordinary samples, no need to enable this mode [experimental, not stable].
 Write additional files to speed up the debugging of graph cleaning.
 
 
+.SS Hi-C-partition options [experimental, not stable]
+
+.TP
+.BI --h1 \ FILEs 
+File names of input Hi-C R1 [r1_1.fq,r1_2.fq,...]
+
+.TP
+.BI --h2 \ FILEs 
+File names of input Hi-C R2 [r2_1.fq,r2_2.fq,...]
+
+
 .SH OUTPUTS
 
 .PP
@@ -347,6 +384,21 @@ phased maternal/haplotype2 contig graph. This graph keeps the phased
 maternal/haplotype2 assembly.
 .RE
 
+.PP
+With Hi-C partition, hifiasm outputs the assembly graphs like trio partition,
+but with additional prefix
+.B [hic].
+In this mode, hifiasm keeps Hi-C alignment results and Hi-C index in two bin 
+files: 
+.B *hic.lk.bin
+and 
+.B *hic.tlb.bin.
+Rerunning hifiasm with different Hi-C reads needs to delete these bin files.
+.RE
+
+
+
+
 .PP
 For each graph, hifiasm also outputs a simplified version without sequences for
 the ease of visualization. Hifiasm keeps corrected reads and overlaps in three