bug fix

I will not create another PR for this...

Assembly.cpp

@@ -10,6 +10,7 @@
 #include "Correct.h"
 #include "htab.h"
 #include "kthread.h"
+#include "rcut.h"
 
 void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, 
 int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct);
@@ -1605,6 +1606,7 @@ void ha_overlap_final(void)
 
 int ha_assemble(void)
 {
+    // debug_mc_g_t(MC_NAME);
 	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)) {

CommandLines.cpp

@@ -23,13 +23,21 @@ static ko_longopt_t long_options[] = {
 	{ "ex-iter",       ko_required_argument, 308 },
     { "purge-cov",     ko_required_argument, 309 },
     { "pri-range",     ko_required_argument, 310 },
-    { "high-het",      ko_no_argument, 311 },
     { "lowQ",          ko_required_argument, 312 },
 	{ "min-hist-cnt",  ko_required_argument, 313 },
     { "h1",      ko_required_argument, 314 },
     { "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 },
+    { "primary",      ko_no_argument, 320 },
+    { "t-occ",      ko_required_argument, 321 },
+    { "seed",      ko_required_argument, 322 },
+    { "n-perturb",      ko_required_argument, 323 },
+    { "f-perturb",      ko_required_argument, 324 },
+    { "n-hap",      ko_required_argument, 325 },
+    { "n-weight",      ko_required_argument, 326 },
 	{ 0, 0, 0 }
 };
 
@@ -45,55 +53,76 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "Usage: hifiasm [options] <in_1.fq> <in_2.fq> <...>\n");
     fprintf(stderr, "Options:\n");
 	fprintf(stderr, "  Input/Output:\n");
-    fprintf(stderr, "    -o STR      prefix of output files [%s]\n", asm_opt->output_file_name);
-    fprintf(stderr, "    -i          ignore saved read correction and overlaps\n");
-    fprintf(stderr, "    -t INT      number of threads [%d]\n", asm_opt->thread_num);
-    fprintf(stderr, "    -z INT      length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
-    fprintf(stderr, "    --version   show version number\n");
+    fprintf(stderr, "    -o STR       prefix of output files [%s]\n", asm_opt->output_file_name);
+    fprintf(stderr, "    -i           ignore saved read correction and overlaps\n");
+    fprintf(stderr, "    -t INT       number of threads [%d]\n", asm_opt->thread_num);
+    fprintf(stderr, "    -z INT       length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
+    fprintf(stderr, "    --version    show version number\n");
 	fprintf(stderr, "  Overlap/Error correction:\n");
-    fprintf(stderr, "    -k INT      k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
-	fprintf(stderr, "    -w INT      minimizer window size [%d]\n", asm_opt->mz_win);
-	fprintf(stderr, "    -f INT      number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift);
-	fprintf(stderr, "    -D FLOAT    drop k-mers occurring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
-	fprintf(stderr, "    -N INT      consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
-    fprintf(stderr, "    -r INT      round of correction [%d]\n", asm_opt->number_of_round);
+    fprintf(stderr, "    -k INT       k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
+	fprintf(stderr, "    -w INT       minimizer window size [%d]\n", asm_opt->mz_win);
+	fprintf(stderr, "    -f INT       number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift);
+	fprintf(stderr, "    -D FLOAT     drop k-mers occurring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
+	fprintf(stderr, "    -N INT       consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
+    fprintf(stderr, "    -r INT       round of correction [%d]\n", asm_opt->number_of_round);
     fprintf(stderr, "  Assembly:\n");
-    fprintf(stderr, "    -a INT      round of assembly cleaning [%d]\n", asm_opt->clean_round);
-    fprintf(stderr, "    -m INT      pop bubbles of <INT in size in contig graphs [%lld]\n", asm_opt->large_pop_bubble_size);
-    fprintf(stderr, "    -p INT      pop bubbles of <INT in size in unitig graphs [%lld]\n", asm_opt->small_pop_bubble_size);
-    fprintf(stderr, "    -n INT      remove tip unitigs composed of <=INT reads [%d]\n", asm_opt->max_short_tip);
-    fprintf(stderr, "    -x FLOAT    max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate);
-    fprintf(stderr, "    -y FLOAT    min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate);
-    fprintf(stderr, "    -u          disable post join contigs step which may improve N50\n");
-    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, "    -a INT       round of assembly cleaning [%d]\n", asm_opt->clean_round);
+    fprintf(stderr, "    -m INT       pop bubbles of <INT in size in contig graphs [%lld]\n", asm_opt->large_pop_bubble_size);
+    fprintf(stderr, "    -p INT       pop bubbles of <INT in size in unitig graphs [%lld]\n", asm_opt->small_pop_bubble_size);
+    fprintf(stderr, "    -n INT       remove tip unitigs composed of <=INT reads [%d]\n", asm_opt->max_short_tip);
+    fprintf(stderr, "    -x FLOAT     max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate);
+    fprintf(stderr, "    -y FLOAT     min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate);
+    fprintf(stderr, "    -u           disable post join contigs step which may improve N50\n");
+    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 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, "    --primary    output a primary assembly and an alternate assembly\n");
     
 //	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");
 
     fprintf(stderr, "  Trio-partition:\n");
-    fprintf(stderr, "    -1 FILE     hap1/paternal k-mer dump generated by \"yak count\" []\n");
-    fprintf(stderr, "    -2 FILE     hap2/maternal k-mer dump generated by \"yak count\" []\n");
-    fprintf(stderr, "    -c INT      lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt);
-    fprintf(stderr, "    -d INT      upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt);
-	fprintf(stderr, "    -3 FILE     list of hap1/paternal read names []\n");
-	fprintf(stderr, "    -4 FILE     list of hap2/maternal read names []\n");
+    fprintf(stderr, "    -1 FILE      hap1/paternal k-mer dump generated by \"yak count\" []\n");
+    fprintf(stderr, "    -2 FILE      hap2/maternal k-mer dump generated by \"yak count\" []\n");
+    fprintf(stderr, "    -c INT       lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt);
+    fprintf(stderr, "    -d INT       upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt);
+	fprintf(stderr, "    -3 FILE      list of hap1/paternal read names []\n");
+	fprintf(stderr, "    -4 FILE      list of hap2/maternal read names []\n");
+    fprintf(stderr, "    --t-occ      INT\n");
+    fprintf(stderr, "                 force remove unitigs with >INT unexpected haplotype-specific reads;\n");
+    fprintf(stderr, "                 ignore graph topology; [%d]\n", asm_opt->trio_flag_occ_thres);
 
     fprintf(stderr, "  Purge-dups:\n");
-    fprintf(stderr, "    -l INT      purge level. 0: no purging; 1: light; 2: aggressive [0 for trio; 2 for unzip]\n");
-    fprintf(stderr, "    -s FLOAT    similarity threshold for duplicate haplotigs [%g]\n", 
-                                     asm_opt->purge_simi_rate);
-    fprintf(stderr, "    -O INT      min number of overlapped reads for duplicate haplotigs [%d]\n", 
-                                     asm_opt->purge_overlap_len);
-    fprintf(stderr, "    --purge-cov INT\n");
-    fprintf(stderr, "                coverage upper bound of Purge-dups [auto]\n");
-    fprintf(stderr, "    --high-het  enable this mode for high heterozygosity sample [experimental, not stable]\n");
+    fprintf(stderr, "    -l INT       purge level. 0: no purging; 1: light; 2/3: aggressive [0 for trio; 3 for unzip]\n");
+    fprintf(stderr, "    -s FLOAT     similarity threshold for duplicate haplotigs [%g for -l1/-l2, %g for -l3]\n", 
+                                      asm_opt->purge_simi_rate_l2, asm_opt->purge_simi_rate_l3);
+    fprintf(stderr, "    -O INT       min number of overlapped reads for duplicate haplotigs [%d]\n", 
+                                      asm_opt->purge_overlap_len);
+    fprintf(stderr, "    --purge-cov  INT\n");
+    fprintf(stderr, "                 coverage upper bound of Purge-dups [auto]\n");
+    fprintf(stderr, "    --n-hap      INT\n");
+    fprintf(stderr, "                 number of haplotypes [%d]\n", asm_opt->polyploidy);
     
-    fprintf(stderr, "  Hi-C-partition [experimental, not stable]:\n");
+    // fprintf(stderr, "  Hi-C-partition [experimental, not stable]:\n");
+    fprintf(stderr, "  Hi-C-partition:\n");
     fprintf(stderr, "    --h1 FILEs   file names of Hi-C R1  [r1_1.fq,r1_2.fq,...]\n");
     fprintf(stderr, "    --h2 FILEs   file names of Hi-C R2  [r2_1.fq,r2_2.fq,...]\n");
+    fprintf(stderr, "    --seed INT   RNG seed [%lu]\n", asm_opt->seed);
+
+    
+    fprintf(stderr, "    --n-weight   INT\n");
+    fprintf(stderr, "                 rounds of reweighting Hi-C links [%d]\n", asm_opt->n_weight);
+	fprintf(stderr, "    --n-perturb  INT\n");
+    fprintf(stderr, "                 rounds of perturbation [%d]\n", asm_opt->n_perturb);
+    fprintf(stderr, "    --f-perturb  FLOAT\n");
+	fprintf(stderr, "                 fraction to flip for perturbation [%.3g]\n", asm_opt->f_perturb);
+
 
     fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
     fprintf(stderr, "See `man ./hifiasm.1' for detailed description of these command-line options.\n");
@@ -102,7 +131,8 @@ void Print_H(hifiasm_opt_t* asm_opt)
 void init_opt(hifiasm_opt_t* asm_opt)
 {
 	memset(asm_opt, 0, sizeof(hifiasm_opt_t));
-	asm_opt->flag = 0;
+	///asm_opt->flag = 0;
+    asm_opt->flag = HA_F_PARTITION;
     asm_opt->coverage = -1;
     asm_opt->num_reads = 0;
     asm_opt->read_file_names = NULL;
@@ -128,7 +158,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->number_of_round = 3;
     asm_opt->adapterLen = 0;
     asm_opt->clean_round = 4;
-    asm_opt->small_pop_bubble_size = 100000;
+    ///asm_opt->small_pop_bubble_size = 100000;
+    asm_opt->small_pop_bubble_size = 0;
     asm_opt->large_pop_bubble_size = 10000000;
     asm_opt->min_drop_rate = 0.2;
     asm_opt->max_drop_rate = 0.8;
@@ -140,20 +171,31 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->max_short_tip = 3;
     asm_opt->min_cnt = 2;
     asm_opt->mid_cnt = 5;
-    asm_opt->purge_level_primary = 2;
+    asm_opt->purge_level_primary = 3;
     asm_opt->purge_level_trio = 0;
-    asm_opt->purge_simi_rate = 0.75;
-    asm_opt->purge_simi_rate_hic = 0.85;
+    asm_opt->purge_simi_rate_l2 = 0.75;
+    asm_opt->purge_simi_rate_l3 = 0.55;
     asm_opt->purge_overlap_len = 1;
-    asm_opt->purge_overlap_len_hic = 50;
+    ///asm_opt->purge_overlap_len_hic = 50;
     asm_opt->recover_atg_cov_min = -1024;
     asm_opt->recover_atg_cov_max = INT_MAX;
     asm_opt->hom_global_coverage = -1;
+    asm_opt->hom_global_coverage_set = 0;
     asm_opt->bed_inconsist_rate = 70;
     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;
+    asm_opt->hap_occ = 1;
+    asm_opt->polyploidy = 2;
+    asm_opt->trio_flag_occ_thres = 60;
+    asm_opt->seed = 11;
+    asm_opt->n_perturb = 10000;
+    asm_opt->f_perturb = 0.1;
+    asm_opt->n_weight = 3;
+    asm_opt->is_alt = 0;
 }
 
 void destory_enzyme(enzyme* f)
@@ -350,9 +392,9 @@ int check_option(hifiasm_opt_t* asm_opt)
         return 0;
     }
 
-    if(asm_opt->purge_level_primary < 0 || asm_opt->purge_level_primary > 2)
+    if(asm_opt->purge_level_primary < 0 || asm_opt->purge_level_primary > 3)
     {
-        fprintf(stderr, "[ERROR] the level of purge-dup should be [0, 2] (-l)\n");
+        fprintf(stderr, "[ERROR] the level of purge-dup should be [0, 3] (-l)\n");
         return 0;
     }
 
@@ -419,25 +461,35 @@ int check_option(hifiasm_opt_t* asm_opt)
         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);
-    // fprintf(stderr, "number of rounds for correction: %d\n", asm_opt->number_of_round);
-    // fprintf(stderr, "number of rounds for assembly cleaning: %d\n", asm_opt->clean_round);
-    // fprintf(stderr, "length of removed adapters: %d\n", asm_opt->adapterLen);
-    // fprintf(stderr, "length of k_mer: %d\n", asm_opt->k_mer_length);
-    // fprintf(stderr, "min overlap drop ratio: %.2g\n", asm_opt->min_drop_rate);
-    // fprintf(stderr, "max overlap drop ratio: %.2g\n", asm_opt->max_drop_rate);
-    // fprintf(stderr, "size of popped small bubbles: %lld\n", asm_opt->small_pop_bubble_size);
-    // fprintf(stderr, "size of popped large bubbles: %lld\n", asm_opt->large_pop_bubble_size);
-    // fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip);
-    // fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip);
-    // fprintf(stderr, "min_cnt: %d\n", asm_opt->min_cnt);
-    // fprintf(stderr, "mid_cnt: %d\n", asm_opt->mid_cnt);
-    // fprintf(stderr, "purge_level_primary: %d\n", asm_opt->purge_level_primary);
-    // fprintf(stderr, "purge_level_trio: %d\n", asm_opt->purge_level_trio);
-    // fprintf(stderr, "purge_simi_rate: %f\n", asm_opt->purge_simi_rate);
-    // fprintf(stderr, "purge_overlap_len: %d\n", asm_opt->purge_overlap_len);
+    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;
+    }
+
+    if(asm_opt->purge_simi_thres < 0)
+    {
+        fprintf(stderr, "[ERROR] [-s] must >= 0\n");
+        return 0;
+    }
 
     return 1;
 }
@@ -575,7 +627,11 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 		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 == 309) asm_opt->hom_global_coverage = atoi(opt.arg);
+        else if (c == 309)
+        {
+            asm_opt->hom_global_coverage = atoi(opt.arg);
+            asm_opt->hom_global_coverage_set = 1;
+        } 
         else if (c == 310)
         {
             char* s = NULL;
@@ -586,18 +642,27 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
                 asm_opt->recover_atg_cov_min = asm_opt->recover_atg_cov_max = -1;
             }
         }
-        else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET;
+        ///else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET;
         else if (c == 312) asm_opt->bed_inconsist_rate = atoi(opt.arg);
         else if (c == 313) asm_opt->min_hist_kmer_cnt = atoi(opt.arg);
         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 == 320) asm_opt->flag -= HA_F_PARTITION, asm_opt->is_alt = 1;
+        else if (c == 321) asm_opt->trio_flag_occ_thres = atoi(opt.arg);
+        else if (c == 322) asm_opt->seed = atol(opt.arg);
+        else if (c == 323) asm_opt->n_perturb = atoi(opt.arg);
+        else if (c == 324) asm_opt->f_perturb = atof(opt.arg);
+        else if (c == 325) asm_opt->polyploidy = atoi(opt.arg);
+        else if (c == 326) asm_opt->n_weight = 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);
         }
-        else if (c == 's') asm_opt->purge_simi_rate = atof(opt.arg);
+        else if (c == 's') asm_opt->purge_simi_rate_l2 = asm_opt->purge_simi_rate_l3 = atof(opt.arg);
         else if (c == 'O') asm_opt->purge_overlap_len = atoll(opt.arg);
         else if (c == ':') 
         {
@@ -611,6 +676,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 		}
     }
 
+    if(asm_opt->purge_level_primary > 2) asm_opt->purge_simi_thres = asm_opt->purge_simi_rate_l3;
+    else asm_opt->purge_simi_thres = asm_opt->purge_simi_rate_l2;
+
     
     if (argc == opt.ind)
     {
@@ -621,6 +689,5 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
     get_queries(argc, argv, &opt, asm_opt);
 
 
-
     return check_option(asm_opt);
 }

CommandLines.h

@@ -2,8 +2,9 @@
 #define __COMMAND_LINE_PARSER__
 
 #include <pthread.h>
+#include <stdint.h>
 
-#define HA_VERSION "0.14-r309"
+#define HA_VERSION "0.15.1-r334"
 
 #define VERBOSE 0
 
@@ -18,6 +19,7 @@
 #define HA_F_BAN_POST_JOIN   0x100
 #define HA_F_BAN_ASSEMBLY    0x200
 #define HA_F_HIGH_HET        0x400
+#define HA_F_PARTITION       0x800
 
 #define HA_MIN_OV_DIFF       0.02 // min sequence divergence in an overlap
 
@@ -49,7 +51,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;
@@ -68,18 +72,22 @@ typedef struct {
     int purge_level_primary;
     int purge_level_trio;
     int purge_overlap_len;
-    int purge_overlap_len_hic;
+    ///int purge_overlap_len_hic;
     int recover_atg_cov_min;
     int recover_atg_cov_max;
     int hom_global_coverage;
+    int hom_global_coverage_set;
     int bed_inconsist_rate;
     int hic_inconsist_rate;
 
     float max_hang_rate;
     float min_drop_rate;
     float max_drop_rate;
-    float purge_simi_rate;
-    float purge_simi_rate_hic;
+    float purge_simi_rate_l2;
+    float purge_simi_rate_l3;
+    float purge_simi_thres;
+
+    ///float purge_simi_rate_hic;
 
     long long small_pop_bubble_size;
     long long large_pop_bubble_size;
@@ -88,7 +96,14 @@ typedef struct {
     long long num_recorrected_bases;
 	long long mem_buf;
     long long coverage;
-    
+    int hap_occ;
+    int polyploidy;
+    int trio_flag_occ_thres;
+    uint64_t seed;
+    int32_t n_perturb;
+    double f_perturb;
+    int32_t n_weight;
+    uint32_t is_alt;
 } hifiasm_opt_t;
 
 extern hifiasm_opt_t asm_opt;

Makefile

@@ -6,7 +6,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 extract.o sys.o ksw2_extz2_sse.o hic.o
+			htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o rcut.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm
 
@@ -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
@@ -72,3 +72,4 @@ main.o: Levenshtein_distance.h htab.h
 sketch.o: kvec.h htab.h Process_Read.h Overlaps.h kdq.h CommandLines.h
 sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
 hic.o: hic.h
+rcut.o: rcut.h

Overlaps.h

@@ -52,6 +52,8 @@
 #define CUT_DIF_HAP 12
 
 
+
+
 ///query is the read itself
 typedef struct {
 	uint64_t qns;
@@ -98,6 +100,35 @@ int max_hang, int min_ovlp);
 long long get_specific_overlap(ma_hit_t_alloc* x, uint32_t qn, uint32_t tn);
 
 
+typedef struct {
+	uint32_t qSpre, qEpre, qScur, qEcur, qn;///[qSp, qEp) && [qSn, qEn]
+	uint32_t tSpre, tEpre, tScur, tEcur, tn;
+} u_trans_hit_t;
+
+typedef struct {
+	size_t n, m;
+	u_trans_hit_t* a;
+} kv_u_trans_hit_t;
+
+
+
+typedef struct {
+	uint32_t qs, qe, qn;
+	uint32_t ts, te, tn;
+	uint32_t occ;
+	double nw;
+	uint8_t f:6, rev:1, del:1;
+} u_trans_t;
+
+typedef struct {
+	size_t n, m;
+	u_trans_t* a;
+	kvec_t(uint64_t) idx;
+} kv_u_trans_t;
+
+#define u_trans_a(x, id) ((x).a + ((x).idx.a[(id)]>>32))
+#define u_trans_n(x, id) ((uint32_t)((x).idx.a[(id)]))
+
 typedef struct {
 	uint64_t ul;
 	uint32_t v;
@@ -289,10 +320,11 @@ static inline asg_arc_t *asg_arc_pushp(asg_t *g)
 // set asg_arc_t::del for v->w
 static inline void asg_arc_del(asg_t *g, uint32_t v, uint32_t w, int del)
 {
-	uint32_t i, nv = asg_arc_n(g, v);
+	uint32_t i, nv = asg_arc_n(g, v)/**, found = 0**/;
 	asg_arc_t *av = asg_arc_a(g, v);
 	for (i = 0; i < nv; ++i)
-		if (av[i].v == w) av[i].del = !!del;
+		if (av[i].v == w) av[i].del = !!del/**, found = 1**/;
+	/**if(found == 0) fprintf(stderr, "ERROR\n");**/
 }
 
 // set asg_arc_t::del and asg_seq_t::del to 1 for sequence s and all its associated arcs
@@ -470,6 +502,7 @@ uint64_t* source_index, long long listLen);
 typedef struct {
 	uint64_t len;
 	uint32_t* index;
+	uint8_t* is_het;
 } R_to_U;
 
 void init_R_to_U(R_to_U* x, uint64_t len);
@@ -477,8 +510,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);
-int asg_pop_bubble_primary(asg_t *g, int max_dist);
+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);
 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);
 
 typedef struct {
@@ -493,9 +525,21 @@ typedef struct {
 	uint32_t new_edges_i;
 } Edge_iter;
 
+typedef struct {
+	asg_arc_t x;
+    uint64_t Off;
+    uint64_t weight;
+}asg_arc_t_offset;
+
+typedef struct {
+	kvec_t(asg_arc_t_offset) a;
+	uint64_t i;
+}kvec_asg_arc_t_offset;
+
+
+
 void init_Edge_iter(asg_t* g, uint32_t v, asg_arc_t* new_edges, uint32_t new_edges_n, Edge_iter* x);
 int get_arc_t(Edge_iter* x, asg_arc_t* get);
-int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag);
 
 
 inline int get_real_length(asg_t *g, uint32_t v, uint32_t* v_s)
@@ -540,55 +584,6 @@ inline uint32_t check_tip(asg_t *sg, uint32_t begNode, uint32_t* endNode, buf_t*
     }
 }
 
-inline uint32_t get_unitig_back(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode, 
-long long* nodeLen, long long* baseLen, buf_t* b)
-{
-	ma_utg_v* u = NULL; 
-    uint32_t v = begNode, w, k;
-    uint32_t kv;
-	(*nodeLen) = (*baseLen) = 0;
-    (*endNode) = (uint32_t)-1;
-	if(ug!=NULL) u = &(ug->u);
-
-    while (1)
-    {
-        kv = get_real_length(sg, v, NULL);
-        (*endNode) = v;
-		if(u == NULL) 
-		{
-			(*nodeLen)++;
-		}
-		else
-		{
-			(*nodeLen) += EvaluateLen((*u), v>>1);
-		}
-        if(b) kv_push(uint32_t, b->b, v);
-		///means reach the end of a unitig
-		if(kv!=1) (*baseLen) += sg->seq[v>>1].len;
-		if(kv==0) return END_TIPS;
-        if(kv>1) return MUL_OUTPUT;
-        ///kv must be 1 here
-        kv = get_real_length(sg, v, &w);
-		///means reach the end of a unitig
-        if(get_real_length(sg, w^1, NULL)!=1)
-		{
-			(*baseLen) += sg->seq[v>>1].len;
-			return MUL_INPUT;
-		}
-
-		for (k = 0; k < asg_arc_n(sg, v); k++)
-        {
-			if(asg_arc_a(sg, v)[k].del) continue;
-			///here is just one undeleted edge
-			(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
-			break;
-		}
-
-        v = w;
-        if(v == begNode) return LOOP;
-    }
-}
-
 inline uint32_t get_unitig(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode, 
 long long* nodeLen, long long* baseLen, long long* max_stop_nodeLen, long long* max_stop_baseLen, 
 uint32_t stops_threshold, buf_t* b)
@@ -1019,36 +1014,55 @@ typedef struct {
 	uint32_t total;
 } Trio_counter;
 
+typedef struct {
+	uint32_t p; // the optimal parent vertex
+	uint32_t d; // the shortest distance from the initial vertex
+	uint32_t r:31, s:1; // r: the number of remaining incoming arc; s: state
+} binfo_s_t;
+
+typedef struct {
+	///all information for each node
+	binfo_s_t *a;
+	kvec_t(uint32_t) S; // set of vertices without parents, nodes with all incoming edges visited
+	kvec_t(uint32_t) b; // visited vertices
+	kvec_t(uint32_t) e; // visited edges/arcs
+} buf_s_t;
+
+typedef struct{
+    buf_s_t *b;
+    uint32_t n_thres, n_reads;
+    asg_t *g;
+    uint32_t check_cross;
+	uint64_t bub_dist;
+} bub_label_t;
+
 void resolve_tangles(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig, 
 long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t trio_flag, 
 float drop_ratio);
-void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
+void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t);
 void rescue_contained_reads_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
-R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t chainLenThres, uint32_t is_bubble_check, 
-uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes);
+R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t chainLenThres, uint32_t is_bubble_check, 
+uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes, bub_label_t* b_mask_t);
 void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
-R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t is_bubble_check, 
-uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges);
+R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t is_bubble_check, uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t);
 void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, 
-ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate);
+ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate, int non_tig_occ);
 void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
 void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g,  ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
-ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges);
+ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges, bub_label_t* b_mask_t);
 void get_unitig_trio_flag(ma_utg_t* nsu, uint32_t flag, uint32_t* require, uint32_t* non_require, uint32_t* ambigious);
 void rescue_missing_overlaps_backward(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
-R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t backward_steps, 
-uint32_t is_bubble_check, uint32_t is_primary_check);
+R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t backward_steps, uint32_t is_bubble_check, uint32_t is_primary_check, bub_label_t* b_mask_t);
 uint32_t get_edge_from_source(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, 
 R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t query, uint32_t target, asg_arc_t* t);
-uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, int max_dist, buf_t *b, 
-uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes);
 int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio);
-
+void asg_bub_backtrack_primary(asg_t *g, uint32_t v0, buf_t *b);
 
 typedef struct{
     double weight;
-    uint32_t uID:31, del:1;
+    uint32_t uID;
     uint64_t dis;
+	uint8_t is_cc:7, del:1;
 	uint64_t occ;
 	///uint64_t occ:63, scaff:1;
     ///uint32_t enzyme;
@@ -1071,35 +1085,87 @@ typedef struct{
 typedef struct{
     kvec_t(hc_linkeage) a;
     kvec_t(uint64_t) enzymes;
-	kvec_t(bed_in) bed;
-    uint32_t* u_idx;
-	uint64_t r_num;
 } hc_links;
 
+#define N_HET 0
+#define C_HET 1
+#define P_HET 2
+#define S_HET 4
+
+typedef struct {
+	uint32_t p_x_p, p_y_p, p_x, p_y;
+	uint32_t c_x_p, c_y_p;
+	uint8_t c_rev;
+} ca_buf_t;
+
+typedef struct {
+	size_t n, m;
+	ca_buf_t* a;
+} kv_ca_buf_t;
+
+typedef struct {
+	kvec_t(uint32_t) uIDs;
+	kvec_t(uint32_t) iDXs;
+	uint32_t chain_num;
+} sub_tran_t;
+
+typedef struct{
+	uint32_t* rUidx;
+	uint64_t* rUpos;
+	uint8_t* is_r_het;
+	uint32_t r_num, u_num;
+	kvec_t(bed_in) bed;
+	kvec_t(uint32_t) topo_buf;
+	kvec_t(uint32_t) topo_res;
+	buf_t b_buf_0, b_buf_1;
+	///uint32_t* uLen;
+	kv_u_trans_t k_trans;
+	kv_u_trans_hit_t k_t_b;
+	kv_ca_buf_t c_buf;
+	sub_tran_t st;
+}trans_chain;
+
+typedef struct {
+    uint32_t n;
+    uint32_t* cov;
+    uint64_t* pos_idx;
+    ma_hit_t_alloc* reverse_sources;
+    ma_sub_t *coverage_cut;
+    R_to_U* ruIndex;
+	asg_t *read_g;
+    int max_hang;
+    int min_ovlp;
+    kvec_asg_arc_t_offset u_buffer; 
+    kvec_t_i32_warp tailIndex;
+    kvec_t_i32_warp prevIndex;
+	///hc_links* link;
+	trans_chain* t_ch;
+}hap_cov_t;
+
 typedef struct{
     ///kvec_t(hc_edge) a;
     size_t n, m; 
     hc_edge *a;
 }hc_edge_warp;
 
-void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num);
+void init_hc_links(hc_links* link, uint64_t ug_num, trans_chain* t_ch);
 void destory_hc_links(hc_links* link);
-void clean_primary_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
-long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
-R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen, 
-uint32_t miniBiGraph, float chimeric_rate, int is_final_clean, int just_bubble_pop, 
-float drop_ratio, hc_links* link);
+uint64_t get_bub_pop_max_dist(asg_t *g, buf_t *b);
+uint64_t get_bub_pop_max_dist_advance(asg_t *g, buf_t *b);
+int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov, uint32_t is_update_chain);
+uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, uint64_t max_dist, buf_t *b, uint32_t positive_flag, 
+uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes, hap_cov_t *cov, uint32_t is_update_chain, uint32_t keep_d);
+
 void adjust_utg_by_primary(ma_ug_t **ug, asg_t* read_g, float drop_rate,
 ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, 
-long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
+long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
 R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
-kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link);
-void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg);
+kvec_asg_arc_t_warp* new_rtg_edges, hap_cov_t **i_cov, bub_label_t* b_mask_t, uint32_t collect_p_trans, uint32_t collect_p_trans_f);
 ma_ug_t* copy_untig_graph(ma_ug_t *src);
 ma_ug_t* output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, 
-uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist, 
+uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
 long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, 
-float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench);
+float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench, bub_label_t* b_mask_t);
 asg_t* copy_read_graph(asg_t *src);
 ma_ug_t *ma_ug_gen(asg_t *g);
 void ma_ug_destroy(ma_ug_t *ug);
@@ -1112,6 +1178,21 @@ inline int inter_interval(int a_s, int a_e, int b_s, int b_e, int* i_s, int* i_e
     return 1;
 }
 
+inline uint32_t get_origin_uid(uint32_t v, trans_chain* t_ch, uint32_t *off, uint32_t *idx)
+{
+	if(off) (*off) = (t_ch->rUpos[v>>1]>>32);
+	if(idx) (*idx) = (uint32_t)(t_ch->rUpos[v>>1]);
+    if(t_ch->rUpos[v>>1] == (uint64_t)-1) return (uint32_t)-1;
+    return (uint32_t)(((t_ch->rUidx[v>>1]>>1)<<1) + ((t_ch->rUidx[v>>1]^v)&1));
+}
+void chain_origin_trans_uid_by_distance(hap_cov_t *cov, asg_t *read_sg, 
+uint32_t *pri_a, uint32_t pri_n, uint32_t pri_beg, uint64_t *i_pri_len, 
+uint32_t *aux_a, uint32_t aux_n, uint32_t aux_beg, uint64_t *i_aux_len,
+ma_ug_t *ug, uint32_t flag, double overall_score, const char* cmd);
+int asg_arc_del_trans(asg_t *g, int fuzz);
+void kt_u_trans_t_idx(kv_u_trans_t *ta, uint32_t n);
+uint32_t get_u_trans_spec(kv_u_trans_t *ta, uint32_t qn, uint32_t tn, u_trans_t **r_a, uint32_t *occ);
+
 #define JUNK_COV 5
 #define DISCARD_RATE 0.8
 

Process_Read.h

@@ -101,6 +101,9 @@ typedef struct
 #define MIX_TRIO 3
 #define NON_TRIO 4
 #define DROP 5
+#define SET_TRIO 8
+#define CHAIN_MATCH 1
+#define CHAIN_UNMATCH 0.334
 
 typedef struct
 {

Purge_Dups.h

@@ -11,14 +11,83 @@
 #define HET_PEAK_RATE (HOM_PEAK_RATE*2)
 #define ALTER_COV_THRES 0.9
 #define REAL_ALTER_THRES 0.1
+#define CHAIN_FILTER_RATE 0.7
+
+#define SELF_EXIST 0
+#define REVE_EXIST 1
+#define DELETE 2
+#define MIXED 3
+#define FLIP 4
+
+#define X2Y 0
+#define Y2X 1
+#define XCY 2
+#define YCX 3
+
+#define Cal_Off(OFF) ((long long)((uint32_t)((OFF)>>32)) - (long long)((uint32_t)((OFF))))
+#define Get_xOff(OFF) ((long long)((uint32_t)((OFF)>>32)))
+#define Get_yOff(OFF) ((long long)((uint32_t)((OFF))))
+#define Get_match(x) ((x).weight)
+#define Get_total(x) ((x).index_beg)
+#define Get_type(x) ((x).index_end)
+#define Get_x_beg(x) ((x).x_beg_pos)
+#define Get_x_end(x) ((x).x_end_pos)
+#define Get_y_beg(x) ((x).y_beg_pos)
+#define Get_y_end(x) ((x).y_end_pos)
+#define Get_rev(x) ((x).rev)
+
+typedef struct {
+    uint8_t rev;
+    uint8_t type;
+    uint8_t status;
+    uint32_t x_beg_pos;
+    uint32_t x_end_pos;
+    uint32_t y_beg_pos;
+    uint32_t y_end_pos;
+    uint32_t x_beg_id;
+    uint32_t x_end_id;
+    uint32_t y_beg_id;
+    uint32_t y_end_id;
+    uint32_t xUid;
+    uint32_t yUid;
+    uint32_t weight;
+    long long score;
+    float s;
+}hap_overlaps;
+
+typedef struct {
+    kvec_t(hap_overlaps) a;
+}kvec_hap_overlaps;
+
+typedef struct {
+    kvec_hap_overlaps* x;
+    uint32_t num;
+}hap_overlaps_list;
 
 void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, 
 ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, 
-uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, 
-uint32_t just_contain, uint32_t just_coverage, hc_links* link);
+uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain, 
+uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t collect_p_trans_f);
 void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge,
 uint32_t is_circle, uint64_t* rLen);
 void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen);
 void enable_debug_mode(uint32_t mode);
+hap_cov_t* init_hap_cov_t(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* sources, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, 
+ma_sub_t *coverage_cut, int max_hang, int min_ovlp, uint32_t is_collect_trans);
+void destory_hap_cov_t(hap_cov_t **x);
+void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd);
+int get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn);
+void set_reverse_hap_overlap(hap_overlaps* dest, hap_overlaps* source, uint32_t* types);
+void print_hap_paf(ma_ug_t *ug, hap_overlaps* ovlp);
+uint64_t get_xy_pos_by_pos(asg_t *read_g, asg_arc_t* t, uint32_t v_in_unitig, uint32_t w_in_unitig, 
+uint32_t v_in_pos, uint32_t w_in_pos, uint32_t xUnitigLen, uint32_t yUnitigLen, uint8_t* rev);
+void quick_LIS(asg_arc_t_offset* x, uint32_t n, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex);
+uint32_t classify_hap_overlap(long long xBeg, long long xEnd, long long xLen,
+long long yBeg, long long yEnd, long long yLen, long long* r_xBeg, long long* r_xEnd, 
+long long* r_yBeg, long long* r_yEnd);
+int cmp_hap_alignment_chaining(const void * a, const void * b);
+uint32_t classify_hap_overlap(long long xBeg, long long xEnd, long long xLen,
+long long yBeg, long long yEnd, long long yLen, long long* r_xBeg, long long* r_xEnd, 
+long long* r_yBeg, long long* r_yEnd);
 
 #endif

README.md

@@ -1,4 +1,4 @@
-## Getting Started
+## <a name="started"></a>Getting Started
 
 ```sh
 # Install hifiasm (requiring g++ and zlib)
@@ -12,30 +12,47 @@ awk '/^S/{print ">"$2;print $3}' test.p_ctg.gfa > test.p_ctg.fa  # get primary c
 
 # Assemble inbred/homozygous genomes (-l0 disables duplication purging)
 hifiasm -o CHM13.asm -t32 -l0 CHM13-HiFi.fa.gz 2> CHM13.asm.log
-# Assemble heterozygous with built-in duplication purging
+# Assemble heterozygous genomes with built-in duplication purging
 hifiasm -o HG002.asm -t32 HG002-file1.fq.gz HG002-file2.fq.gz
 
+# Hi-C phasing with paired-end short reads in two FASTQ files
+hifiasm -o HG002.asm --h1 read1.fq.gz --h2 read2.fq.gz HG002-HiFi.fq.gz
+
 # Trio binning assembly (requiring https://github.com/lh3/yak)
 yak count -b37 -t16 -o pat.yak <(cat pat_1.fq.gz pat_2.fq.gz) <(cat pat_1.fq.gz pat_2.fq.gz)
 yak count -b37 -t16 -o mat.yak <(cat mat_1.fq.gz mat_2.fq.gz) <(cat mat_1.fq.gz mat_2.fq.gz)
 hifiasm -o HG002.asm -t32 -1 pat.yak -2 mat.yak HG002-HiFi.fa.gz
 ```
 
-## Introduction
+## Table of Contents
 
-Hifiasm is a fast haplotype-resolved de novo assembler for PacBio Hifi reads. 
-It can assemble a human genome in several hours and works with the California
-redwood genome, one of the most complex genomes sequenced so far. Hifiasm can
-produce primary/alternate assemblies of quality competitive with the best
-assemblers. It also introduces a new graph binning algorithm and achieves
-the best haplotype-resolved assembly given trio data.
+- [Getting Started](#started)
+- [Introduction](#intro)
+- [Why Hifiasm?](#why)
+- [Usage](#use)
+  - [Assembling HiFi reads without additional data types](#hifionly)
+  - [Hi-C integration](#hic)
+  - [Trio binning](#trio)
+  - [Output files](#output)
+- [Results](#results)
+- [Getting Help](#help)
+- [Limitations](#limit)
+- [Citing Hifiasm](#cite)
 
-## Why Hifiasm?
+## <a name="intro"></a>Introduction
+
+Hifiasm is a fast haplotype-resolved de novo assembler for PacBio HiFi reads.
+It can assemble a human genome in several hours and assemble a ~30Gb California
+redwood genome in a few days. Hifiasm emits partially phased assemblies of
+quality competitive with the best assemblers. Given parental short reads or
+Hi-C data, it produces arguably the best haplotype-resolved assemblies so far.
+
+## <a name="why"></a>Why Hifiasm?
 
 * Hifiasm delivers high-quality assemblies. It tends to generate longer contigs
   and resolve more segmental duplications than other assemblers.
 
-* Given sequence reads from the parents, hifiasm can produce overall the best
+* Given Hi-C reads or short reads from the parents, hifiasm can produce overall the best
   haplotype-resolved assembly so far. It is the assembler of choice by the
   [Human Pangenome Project][hpp] for the first batch of samples.
 
@@ -47,13 +64,15 @@ the best haplotype-resolved assembly given trio data.
 * Hifiasm is fast. It can assemble a human genome in half a day and assemble a
   ~30Gb redwood genome in three days. No genome is too large for hifiasm.
 
-* Hifiasm is trivial to install and easy to use. It does not required python,
-  R or C++11 compilers and can be compiled into a single executable. The
+* Hifiasm is trivial to install and easy to use. It does not required Python,
+  R or C++11 compilers, and can be compiled into a single executable. The
   default setting works well with a variety of genomes.
 
 [hpp]: https://humanpangenome.org
 
-## Usage
+## <a name="use"></a>Usage
+
+### <a name="hifionly"></a>Assembling HiFi reads without additional data types
 
 A typical hifiasm command line looks like:
 ```sh
@@ -61,11 +80,21 @@ hifiasm -o NA12878.asm -t 32 NA12878.fq.gz
 ```
 where `NA12878.fq.gz` provides the input reads, `-t` sets the number of CPUs in
 use and `-o` specifies the prefix of output files. For this example, the
-primary contigs are written to `NA12878.asm.p_ctg.gfa` and alternate contigs to
-`NA12878.asm.a_ctg.gfa`. At the first run, hifiasm saves corrected reads and
+primary contigs are written to `NA12878.asm.bp.p_ctg.gfa` and alternate contigs to
+`NA12878.asm.bp.a_ctg.gfa`. Since v0.15, hifiasm also produces two sets of
+partially phased contigs at `NA12878.asm.bp.hap?.p_ctg.gfa`. This pair of files
+can be thought to represent the two haplotypes in a diploid genome, though with
+occasional switch errors. The frequency of switches is determined by the
+heterozygosity of the input sample.
+
+At the first run, hifiasm saves corrected reads and
 overlaps to disk as `NA12878.asm.*.bin`. It reuses the saved results to avoid
 the time-consuming all-vs-all overlap calculation next time. You may specify
 `-i` to ignore precomputed overlaps and redo overlapping from raw reads.
+You can also dump error corrected reads in FASTA and read overlaps in PAF with
+```sh
+hifiasm -o NA12878.asm -t 32 --write-paf --write-ec /dev/null
+```
 
 Hifiasm purges haplotig duplications by default. For inbred or homozygous
 genomes, you may disable purging with option `-l0`. Old HiFi reads may contain
@@ -75,7 +104,27 @@ bloom filter which takes 16GB memory at the beginning. For genomes much larger
 than human, applying `-f38` or even `-f39` is preferred to save memory on k-mer
 counting.
 
-When parental short reads are available, hifiasm can generate a pair of
+### <a name="hic"></a>Hi-C integration
+
+Hifiasm can generate a pair of haplotype-resolved assemblies with paired-end
+Hi-C reads:
+```sh
+hifiasm -o NA12878.asm -t32 --h1 read1.fq.gz --h2 read2.fq.gz HiFi-reads.fq.gz
+```
+In this mode, each contig is supposed to be a haplotig, which by definition
+comes from one parental haplotype only. Hifiasm often puts all contigs from the
+same parental chromosome in one assembly. It has cleanly separated chrX and
+chrY for a human male dataset. Nonetheless, phasing across centromeres is
+challenging. Users should not expect hifiasm to phase entire chromosomes at the
+moment. Also, contigs from different parental chromosomes are randomly mixed as
+it is just not possible to phase across chromosomes with Hi-C.
+
+Hifiasm does not perform scaffolding for now. You need to run a standalone
+scaffolder such as SALSA or 3D-DNA to scaffold phased haplotigs.
+
+### <a name="trio"></a>Trio binning
+
+When parental short reads are available, hifiasm can also generate a pair of
 haplotype-resolved assemblies with trio binning. To perform such assembly, you
 need to count k-mers first with [yak][yak] first and then do assembly:
 ```sh
@@ -85,19 +134,15 @@ hifiasm -o NA12878.asm -t 32 -1 pat.yak -2 mat.yak NA12878.fq.gz
 ```
 Here `NA12878.asm.hap1.p_ctg.gfa` and `NA12878.asm.hap2.p_ctg.gfa` give the two
 haplotype assemblies. In the binning mode, hifiasm does not purge haplotig
-duplications by default. Because hifiasm reuses saved overlaps, you can
+duplicates by default. Because hifiasm reuses saved overlaps, you can
 generate both primary/alternate assemblies and trio binning assemblies with
 ```sh
 hifiasm -o NA12878.asm -t 32 NA12878.fq.gz 2> NA12878.asm.pri.log
 hifiasm -o NA12878.asm -t 32 -1 pat.yak -2 mat.yak /dev/null 2> NA12878.asm.trio.log
 ```
-The second command line will run much faster than the first. You can also dump
-error corrected in FASTA and/or overlaps in PAF with
-```sh
-hifiasm -o NA12878.asm -t 32 --write-paf --write-ec /dev/null
-```
+The second command line will run much faster than the first.
 
-## Output files
+### <a name="output"></a>Output files
 
 For non-trio assembly, hifiasm generates the following files:
 
@@ -126,9 +171,9 @@ For trio assembly, hifiasm generates the following files:
 Hifiasm writes error corrected reads to the *prefix*.ec.bin binary file and
 writes overlaps to *prefix*.ovlp.source.bin and *prefix*.ovlp.reverse.bin.
 
-## Results
+## <a name="results"></a>Results
 
-The following table shows the statistics of several hifiasm primary assemblies:
+The following table shows the statistics of several hifiasm primary assemblies assembled with v0.12:
 
 |<sub>Dataset<sub>|<sub>Size<sub>|<sub>Cov.<sub>|<sub>Asm options<sub>|<sub>CPU time<sub>|<sub>Wall time<sub>|<sub>RAM<sub>|<sub> N50<sub>|
 |:---------------|-----:|-----:|:---------------------|-------:|--------:|----:|----------------:|
@@ -155,7 +200,10 @@ redwood genome in a few days on a single machine. For trio binning assembly:
 |:---------------|-----:|-------:|--------:|----:|----------------:|
 |<sub>[HG00733][HG00733-data], [\[father\]][HG00731-data], [\[mother\]][HG00732-data]</sub>|<sub>&times;33</sub>|<sub>269.1h</sub>|<sub>6.9h</sub>|<sub>135G</sub>|<sub>35.1Mb (paternal), 34.9Mb (maternal)</sub>|
 |<sub>[HG002][NA24385-data],   [\[father\]][NA24149-data], [\[mother\]][NA24143-data]</sup>|<sub>&times;36</sub>|<sub>305.4h</sub>|<sub>7.7h</sub>|<sub>137G</sub>|<sub>41.0Mb (paternal), 40.8Mb (maternal)</sub>|
+
+<!--
 |<sub>[NA12878][NA12878-data], [\[father\]][NA12891-data], [\[mother\]][NA12892-data]</sub>|<sub>&times;30</sub>|<sub>180.8h</sub>|<sub>4.9h</sub>|<sub>123G</sub>|<sub>27.7Mb (paternal), 27.0Mb (maternal)</sub>|
+-->
 
 [HG00733-data]: https://www.ebi.ac.uk/ena/data/view/ERX3831682
 [HG00731-data]: https://www.ebi.ac.uk/ena/data/view/ERR3241754
@@ -167,29 +215,32 @@ redwood genome in a few days on a single machine. For trio binning assembly:
 [NA12891-data]: https://www.ebi.ac.uk/ena/data/view/ERR194160
 [NA12892-data]: https://www.ebi.ac.uk/ena/data/view/ERR194161
 
-Except NA12878, the assemblies above were produced by hifiasm v0.12 and can be
-downloaded at
-```txt
-ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/submission/hifiasm-0.12/
-```
-NA12878 was assembled with an older version of hifiasm and is available at
-```txt
-ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/NA12878-r253/
-```
-
+Human assemblies above can be acquired [from Zenodo][zenodo-human] and
+non-human ones are available [here][zenodo-nonh].
 
+[zenodo-human]: https://zenodo.org/record/4393631
+[zenodo-nonh]: https://zenodo.org/record/4393750
 [unitig]: http://wgs-assembler.sourceforge.net/wiki/index.php/Celera_Assembler_Terminology
 [gfa]: https://github.com/pmelsted/GFA-spec/blob/master/GFA-spec.md
 [paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
 [yak]: https://github.com/lh3/yak
 
-## Getting Help
+## <a name="help"></a>Getting Help
 
 For detailed description of options, please see `man ./hifiasm.1`. The `-h`
 option of hifiasm also provides brief description of options. If you have
 further questions, please raise an issue at the [issue
 page](https://github.com/chhylp123/hifiasm/issues).
 
-## Limitations
+## <a name="limit"></a>Limitations
 
-1. Purging haplotig duplications may introduce misassemblies.
+1. Purging haplotig duplications may introduce misassemblies.
+
+## <a name="cite"></a>Citating Hifiasm
+
+If you use hifiasm in your work, please cite:
+
+> Cheng, H., Concepcion, G.T., Feng, X., Zhang, H., Li H. (2021)
+> Haplotype-resolved de novo assembly using phased assembly graphs with
+> hifiasm. *Nat Methods*, **18**:170-175.
+> https://doi.org/10.1038/s41592-020-01056-5

hic.h

@@ -10,8 +10,8 @@
 #define RC_2 2
 
 hc_edge* get_hc_edge(hc_links* link, uint64_t src, uint64_t dest, uint64_t dir);
-void push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t* d);
-void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link);
+hc_edge* push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t* d);
+void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch);
 void hic_benchmark(ma_ug_t *ug, asg_t* read_g);
 
 typedef struct {
@@ -49,6 +49,11 @@ typedef struct {
     kvec_t(chain_w_type) chain_weight;
     chain_hic_warp c_w;
 } bubble_type;
+
+typedef struct {
+    int8_t *s;
+    uint64_t xs;
+} ps_t;
 #define P_het(B) ((B).num.n)
 #define M_het(B) ((B).num.n + 1)
 // #define IF_BUB(ID, B) ((B).index[(ID)] < (B).num.n)
@@ -62,11 +67,12 @@ void get_bubbles(bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t* sink, u
 int load_hc_links(hc_links* link, const char *fn);
 void write_hc_links(hc_links* link, const char *fn);
 void destory_bubbles(bubble_type* bub);
-void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link);
+void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_trans_t *ref);
 void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub);
 uint32_t connect_bub_occ(bubble_type* bub, uint32_t root_id, uint32_t check_het);
 void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, uint32_t check_het);
 void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint32_t is_end);
 void set_b_utg_weight_flag(bubble_type* bub, buf_t* b, uint32_t v, uint8_t* vis_flag, uint32_t flag, uint32_t* occ);
+void debug_gfa_space(ma_ug_t* ug, hap_cov_t *cov);
 
 #endif

hifiasm.1

@@ -1,4 +1,4 @@
-.TH hifiasm 1 "19 July 2020" "hifiasm-0.9 (r289)" "Bioinformatics tools"
+.TH hifiasm 1 "16 April 2021" "hifiasm-0.15 (r327)" "Bioinformatics tools"
 
 .SH NAME
 .PP
@@ -212,6 +212,38 @@ 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.
+
+.TP
+.BI --primary
+Output a primary assembly and an alternate assembly.
+Hifiasm outputs two balanced assemblies and a primary
+assembly in default. Enable this option or 
+.B -l0
+outputs a primary assembly and an alternate assembly.
+
+
 .SS Trio-partition options
 
 .TP 10
@@ -260,12 +292,14 @@ times in the other sample.
 .TP 10
 .BI -l \ INT
 Level of purge-dup. 0 to disable purge-dup, 1 to only purge contained haplotigs, 
-2 to purge all types of haplotigs. In default, [2] for non-trio assembly, [0] for trio assembly.
+2 to purge all types of haplotigs, 3 to purge all types of haplotigs in most aggressive way 
+for high heterozygosity sample. 
+In default, [3] for non-trio assembly, [0] for trio assembly.
 For trio assembly, only level 0 and level 1 are allowed.
 
 .TP
 .BI -s \ FLOAT
-Similarity threshold for duplicate haplotigs that should be purged [0.75].
+Similarity threshold for duplicate haplotigs that should be purged [0.75 for -l1/-l2, 0.55 for -l3].
 
 .TP
 .BI -O \ FLOAT
@@ -277,9 +311,8 @@ Coverage upper bound of Purge-dups, which is inferred automatically in default.
 If the coverage of a contig is higher than this bound, don't apply Purge-dups.
 
 .TP
-.BI --high-het \ INT
-Enable this mode for high heterozygosity sample, which will increase running time. 
-For ordinary samples, no need to enable this mode [experimental, not stable].
+.BI   --n-hap \ INT
+Assumption of haplotype number.
 
 
 .SS Debugging options
@@ -289,14 +322,39 @@ 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,...].
+
+.TP
+.BI --n-weight \ INT 
+Rounds of reweighting Hi-C links [3]. Increasing this may improves
+phasing results but takes longer time.
+
+.TP
+.BI --n-perturb \ INT 
+Rounds of perturbation [10000]. Increasing this may improves
+phasing results but takes longer time.
+
+.TP
+.BI --f-perturb \ FLOAT 
+Fraction to flip for perturbation [0.1]. Increasing this may improves
+phasing results but takes longer time.
+
+.TP
+.BI --seed \ INT 
+RNG seed [11].
+
 .SH OUTPUTS
 
 .PP
-Without trio partition options
-.B -1
-and
-.BR -2 ,
-hifiasm generates the following assembly graphs in the GFA format:
+In general, hifiasm generates the following assembly graphs in the GFA format:
 
 .RS 2
 .TP 2
@@ -323,30 +381,111 @@ assembly graph of primary contigs. This graph collapses different haplotypes.
 assembly graph of alternate contigs. This graph consists of all assemblies that
 are discarded in primary contig graph.
 
+.TP
+*
+.IR prefix .hap*.p_ctg.gfa:
+phased contig graph. This graph keeps the phased assembly.
+
 .RE
 
 .PP
-With trio partition, hifiasm outputs the following assembly graphs:
+Hifiasm outputs
+.B *.r_utg.gfa
+and
+.B *.p_utg.gfa
+in any cases.
+Specifically, hifiasm outputs the following assembly graphs
+with trio-binning options:
 
 .RS 2
 .TP 2
 *
-.IR prefix .dip.r_utg.gfa:
-haplotype-resolved raw unitig graph. This graph keeps all haplotype information.
-
-.TP
-*
-.IR prefix .hap1.p_ctg.gfa:
-phased paternal/haplotype1 contig graph. This graph keeps the phased
+.IR prefix .dip.hap1.p_ctg.gfa:
+phased paternal/haplotype1 contig graph keeping the phased
 paternal/haplotype1 assembly.
 
 .TP
 *
-.IR prefix .hap2.p_ctg.gfa:
-phased maternal/haplotype2 contig graph. This graph keeps the phased
+.IR prefix .dip.hap2.p_ctg.gfa:
+phased maternal/haplotype2 contig graph keeping the phased
 maternal/haplotype2 assembly.
 .RE
 
+.PP
+With Hi-C partition options, hifiasm outputs:
+
+.RS 2
+.TP 2
+*
+.IR prefix .hic.p_ctg.gfa:
+assembly graph of primary contigs. This graph collapses different haplotypes.
+
+.TP
+*
+.IR prefix .hic.hap1.p_ctg.gfa:
+phased contig graph where each contig is fully phased.
+
+.TP
+*
+.IR prefix .hic.hap2.p_ctg.gfa:
+phased contig graph where each contig is fully phased.
+.RE
+
+
+.PP
+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
+or enable 
+.BR -i .
+.RE
+
+.PP
+Hifiasm generates the following assembly graphs only with HiFi reads:
+
+.RS 2
+.TP 2
+*
+.IR prefix .p_ctg.gfa:
+assembly graph of primary contigs. This graph collapses different haplotypes.
+
+.TP
+*
+.IR prefix .bp.hap1.p_ctg.gfa:
+balanced contig graph where each contig is partially phased.
+
+.TP
+*
+.IR prefix .bp.hap2.p_ctg.gfa:
+balanced contig graph where each contig is partially phased.
+.RE
+
+.PP
+If the option
+.BR -p
+or
+.BR --primary
+is specified, hifiasm outputs:
+
+.RS 2
+.TP 2
+*
+.IR prefix .p_ctg.gfa:
+assembly graph of primary contigs. This graph collapses different haplotypes.
+
+.TP
+*
+.IR prefix .a_ctg.gfa:
+assembly graph of alternate contigs. This graph consists of all assemblies that
+are discarded in primary contig graph.
+.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

main.cpp

@@ -11,7 +11,7 @@ int main(int argc, char *argv[])
 	int i, ret;
 	yak_reset_realtime();
     init_opt(&asm_opt);
-    if (!CommandLine_process(argc, argv, &asm_opt)) return 1;
+    if (!CommandLine_process(argc, argv, &asm_opt)) return 0;
 	ret = ha_assemble();
     destory_opt(&asm_opt);
 	fprintf(stderr, "[M::%s] Version: %s\n", __func__, HA_VERSION);

rcut.h

@@ -0,0 +1,93 @@
+#ifndef __RCUT__
+#define __RCUT__
+#include <stdio.h>
+#include <stdint.h>
+#include "kvec.h"
+#include "Overlaps.h"
+#include "Purge_Dups.h"
+#include "hic.h"
+
+typedef struct {
+    uint32_t bS, bE;
+    uint32_t nS, nE;
+    uint32_t uID;
+    uint8_t hs;
+}mc_interval_t;
+
+#define mc_node_t int8_t
+// #define w_t int64_t
+// #define t_w_t int64_t
+// #define w_cast(x) ((t_w_t)((x) < 0 ? (x) - 0.5 : (x) + 0.5))
+
+#define w_t double
+#define t_w_t double
+#define w_cast(x) ((t_w_t)((x)))
+#define MC_NAME "debug_mc.bin"
+
+typedef struct {
+	uint64_t x; ///(uint64_t)nid1 << 32 | nid2;
+	w_t w;  ///might be negative or positive
+} mc_edge_t;
+
+typedef struct {
+    kvec_t(uint64_t) idx;
+    kvec_t(mc_edge_t) ma;
+    uint64_t* cc;
+    uint32_t n_seq;
+} mc_match_t;
+
+typedef struct {
+    kvec_t(mc_node_t) s;
+    ma_ug_t *ug;
+    asg_t *rg;
+    mc_match_t* e;
+}mc_g_t;
+
+typedef struct {
+    uint32_t a[2], occ[2]; 
+    mc_node_t s[2];
+    t_w_t z[4];
+}mb_node_t;
+
+typedef struct {
+    kvec_t(uint32_t) bid;
+    kvec_t(uint32_t) idx;
+    kvec_t(mb_node_t) u;
+}mb_nodes_t;
+
+typedef struct {
+	uint64_t x; ///(uint64_t)nid1 << 32 | nid2;
+	t_w_t w[4];  ///might be negative or positive
+} mb_edge_t;
+
+typedef struct {
+    kvec_t(uint64_t) idx;
+    kvec_t(mb_edge_t) ma;
+    uint64_t* cc;
+    uint32_t n_seq;
+} mb_match_t;
+
+typedef struct {
+    mb_nodes_t* u; 
+    mb_match_t* e;
+}mb_g_t;
+
+static inline uint64_t kr_splitmix64(uint64_t x)
+{
+	uint64_t z = (x += 0x9E3779B97F4A7C15ULL);
+	z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL;
+	z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL;
+	return z ^ (z >> 31);
+}
+
+static inline double kr_drand_r(uint64_t *x)
+{
+    union { uint64_t i; double d; } u;
+	*x = kr_splitmix64(*x);
+    u.i = 0x3FFULL << 52 | (*x) >> 12;
+    return u.d - 1.0;
+}
+
+void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref);
+void debug_mc_g_t(const char* name);
+#endif