add scaffolding

CommandLines.cpp

@@ -62,6 +62,8 @@ static ko_longopt_t long_options[] = {
     { "path-min",     ko_required_argument, 347},
     { "trio-dual",     ko_no_argument, 348},
     { "ul-cut",     ko_required_argument, 349},
+    { "dual-scaf",     ko_no_argument, 350},
+    { "scaf-gap",     ko_required_argument, 351},
     // { "path-round",     ko_required_argument, 348},
 	{ 0, 0, 0 }
 };
@@ -163,7 +165,7 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "    --l-msjoin   INT\n");
     fprintf(stderr, "                 detect misjoined unitigs of >=INT in size; 0 to disable [%lu]\n", asm_opt->misjoin_len);
 
-    fprintf(stderr, "  Ultra-Long-integration (beta):\n");
+    fprintf(stderr, "  Ultra-Long-integration:\n");
     fprintf(stderr, "    --ul FILEs   file names of Ultra-Long reads [r1.fq,r2.fq,...]\n");
     fprintf(stderr, "    --ul-rate    FLOAT\n");
     fprintf(stderr, "                 error rate of Ultra-Long reads [%.3g]\n", asm_opt->ul_error_rate);
@@ -179,6 +181,11 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "                 filter out <INT UL reads during the UL assembly [%d]\n", asm_opt->ul_min_base);
     // fprintf(stderr, "    --low-het    enable it for genomes with very low het heterozygosity rate (<0.0001%%)\n");
 
+    fprintf(stderr, "  Dual-Scaffolding:\n");
+    fprintf(stderr, "    --dual-scaf  output scaffolding\n");
+    fprintf(stderr, "    --scaf-gap   INT\n");
+    fprintf(stderr, "                 max gap size for scaffolding [%ld]\n", asm_opt->self_scaf_gap_max);
+
     fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
     fprintf(stderr, "See `https://hifiasm.readthedocs.io/en/latest/' or `man ./hifiasm.1' for complete documentation.\n");
 }
@@ -295,6 +302,9 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->trio_cov_het_ovlp = -1;
     asm_opt->ul_min_base = 0;
     asm_opt->self_scaf = 0;
+    asm_opt->self_scaf_min = 250000;
+    asm_opt->self_scaf_reliable_min = 5000000;
+    asm_opt->self_scaf_gap_max = 3000000;
 }
 
 void destory_enzyme(enzyme* f)
@@ -841,6 +851,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 347) asm_opt->min_path_drop_rate = atof(opt.arg);
         else if (c == 348) asm_opt->trio_cov_het_ovlp = 1;
         else if (c == 349) asm_opt->ul_min_base = atol(opt.arg);
+        else if (c == 350) asm_opt->self_scaf = 1;
+        else if (c == 351) asm_opt->self_scaf_gap_max = atol(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

@@ -5,7 +5,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.19.6-r597"
+#define HA_VERSION "0.19.7-r598"
 
 #define VERBOSE 0
 
@@ -148,6 +148,9 @@ typedef struct {
     uint8_t hifi_pst_join, ul_pst_join;
     uint32_t ul_min_base;
     uint8_t self_scaf;
+    uint64_t self_scaf_min;
+    uint64_t self_scaf_reliable_min;
+    int64_t self_scaf_gap_max;
 } hifiasm_opt_t;
 
 extern hifiasm_opt_t asm_opt;

Overlaps.h

@@ -1176,7 +1176,7 @@ void extract_sub_overlaps(uint32_t i_tScur, uint32_t i_tEcur, uint32_t i_tSpre,
 uint32_t tn, kv_u_trans_hit_t* ktb, uint32_t bn);
 void clean_u_trans_t_idx(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g);
 void clean_u_trans_t_idx_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g);
-void clean_u_trans_t_idx_filter_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g);
+void clean_u_trans_t_idx_filter_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double sc_sec_rate, uint64_t uniform_only);
 uint32_t test_dbug(ma_ug_t* ug, FILE* fp);
 void write_dbug(ma_ug_t* ug, FILE* fp);
 int asg_arc_identify_simple_bubbles_multi(asg_t *g, bub_label_t* x, int check_cross);

hic.cpp

@@ -17744,7 +17744,7 @@ void trio_phasing_refine(ma_ug_t *iug, asg_t* sg, kv_u_trans_t *ref, ug_opt_t *o
     uint64_t ug_n0 = ug->g->n_seq, k, i, ul, cis_n, trans_n, w_n, tot_hap, tot_r, frid, mrid, flag; 
     kv_u_trans_t in; memset(&in, 0, sizeof(in)); ma_utg_t *p; u_trans_t *z; int64_t fh, mh;
     bubble_type *bub = gen_bubble_chain(sg, ug, opt, &bf, 0); free(bf);
-    clean_u_trans_t_idx_filter_adv(ref, ug, sg);
+    clean_u_trans_t_idx_filter_adv(ref, ug, sg, 0.95, 0);
     
     ///update ug itself
     ul = FATHER; frid = ug->g->n_seq; asg_seq_set(ug->g, frid, ul, 0); 

inter.h

@@ -128,6 +128,8 @@ uint32_t clean_contain_g(const ug_opt_t *uopt, asg_t *sg, uint32_t push_trans);
 void dedup_contain_g(const ug_opt_t *uopt, asg_t *sg);
 void trans_base_mmhap_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res);
 scaf_res_t *gen_contig_path(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *ctg, ma_ug_t *ref);
-kv_u_trans_t *gen_contig_trans(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *qry, scaf_res_t *qry_sc, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, kv_u_trans_t *ta);
+void gen_contig_trans(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *qry, scaf_res_t *qry_sc, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, kv_u_trans_t *ta, uint32_t qoff, uint32_t toff, bubble_type *bu, kv_u_trans_t *res);
+void gen_contig_self(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *db, scaf_res_t *db_sc, ma_ug_t *gfa, kv_u_trans_t *ta, uint64_t soff, bubble_type *bu, kv_u_trans_t *res);
+void order_contig_trans(kv_u_trans_t *in);
 
 #endif