bub disable

Added citation

CommandLines.cpp

@@ -23,13 +23,20 @@ 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 },
 	{ 0, 0, 0 }
 };
 
@@ -45,55 +52,72 @@ 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-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 +126,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 +153,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 +166,30 @@ 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_simi_rate_hic = 0.85;
     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 = 50000;
+    asm_opt->f_perturb = 0.1;
 }
 
 void destory_enzyme(enzyme* f)
@@ -350,9 +386,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 +455,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 +621,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 +636,26 @@ 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;
+        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 == '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 +669,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)
     {

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-r327"
 
 #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,21 @@ 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 +95,12 @@ 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;
 } 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 partig.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,5 @@ 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
+partig.o: partig.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,31 +1014,49 @@ 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;
@@ -1071,35 +1084,80 @@ 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{
+	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;
+}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);
 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 +1170,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,82 @@
 #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;
+}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);
 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

@@ -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.h

@@ -11,7 +11,7 @@
 
 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);
+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,13 @@ 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);
 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);
+uint32_t get_unitig_het_arb(ma_utg_t* u, uint8_t *r_het_flag, uint32_t m_het_label, uint32_t p_het_label, uint32_t n_het_label);
+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,34 @@ 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-perturb \ INT 
+Rounds of perturbation [50000]. Increasing this improves
+phasing results but takes longer time.
+
+.TP
+.BI --f-perturb \ FLOAT 
+Fraction to flip for perturbation [0.1]. Increasing this 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 +376,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

partig.cpp

@@ -0,0 +1,834 @@
+#define __STDC_LIMIT_MACROS
+#include <stdint.h>
+#include <stdlib.h>
+#include "partig.h"
+#include "Purge_Dups.h"
+#include "Correct.h"
+#include "ksort.h"
+
+#define generic_key(x) (x)
+KRADIX_SORT_INIT(gfa64, uint64_t, generic_key, 8)
+
+#define pt_a(x, id) ((x).ma.a + ((x).idx.a[(id)]>>32))
+#define pt_n(x, id) ((uint32_t)((x).idx.a[(id)]))
+
+typedef struct {
+	int32_t topn;
+	int32_t n_perturb;
+	uint64_t seed;
+	double f_perturb;
+} pt_svopt_t;
+
+typedef struct {
+	///uint32_t m, n, *shuffled;
+	uint32_t *shuffled;
+	uint32_t off, size; // offset in pt_match_t::cc; size of the component
+	uint64_t *buf;
+	kvec_t(uint64_t) a;
+	int8_t *s, *s_tmp;///s is the status (haplotype) of each unitig: for backup
+} solve_aux_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 pt_svopt_init(pt_svopt_t *opt)
+{
+	memset(opt, 0, sizeof(pt_svopt_t));
+	opt->seed = 11;
+	opt->topn = 1<<30;
+	opt->n_perturb = 1000;
+	opt->f_perturb = 0.1;
+}
+
+static void pt_pdist_idx(pt_match_t *ma)
+{
+	uint32_t st, i;
+	kv_resize(uint64_t, ma->idx, ma->n_seq);
+	ma->idx.n = ma->n_seq;
+	memset(ma->idx.a, 0, ma->idx.n*sizeof(uint64_t));
+	for (st = 0, i = 1; i <= ma->ma.n; ++i)
+		if (i == ma->ma.n || ma->ma.a[i].sid[0] != ma->ma.a[st].sid[0])
+			ma->idx.a[ma->ma.a[st].sid[0]] = (uint64_t)st << 32 | (i - st), st = i;
+}
+
+static pt_match1_t *pt_pdist(const pt_match_t *ma, uint32_t sid1, uint32_t sid2)
+{
+	pt_match1_t *o = pt_a(*ma, sid1);
+	uint32_t n = pt_n(*ma, sid1), k;
+	for (k = 0; k < n; ++k)
+		if (o[k].sid[1] == sid2)
+			return &(o[k]);
+	return NULL;
+}
+
+static void normalize_pdist(pt_match1_t *a, pt_match1_t *b)
+{
+	if(a->w >= b->w)
+	{
+		b->sid[0] = a->sid[1];
+		b->sid[1] = a->sid[0];
+		b->w = a->w;
+	}
+	else
+	{
+		a->sid[0] = b->sid[1];
+		a->sid[1] = b->sid[0];
+		a->w = b->w;
+	}
+}
+
+uint32_t pt_pdist_symm(pt_match_t *ma)
+{
+	uint8_t *del = NULL;
+	uint32_t i, k, n = 0;
+	pt_match1_t *t = NULL, *m = NULL;
+	CALLOC(del, ma->ma.n);
+
+	for (i = 0; i < ma->ma.n; ++i) {
+		m = &ma->ma.a[i];
+		if (m->sid[0] == m->sid[1])
+		{
+			del[i] = 1, ++n;///self overlap
+			continue;
+		}
+		t = pt_pdist(ma, m->sid[1], m->sid[0]);
+		if(!t)
+		{
+			del[i] = 1, ++n;///self overlap
+			continue;
+		} 
+		normalize_pdist(m, t);
+	}
+
+	if (n > 0) {
+		for (i = k = 0; i < ma->ma.n; ++i)
+			if (!del[i]) ma->ma.a[k++] = ma->ma.a[i];
+		ma->ma.n = k;
+		pt_pdist_idx(ma);
+	}
+
+	free(del);
+	return n;
+}
+
+
+static void pt_pdist_flt(pt_match_t *ma, uint32_t min_cnt, double drop_thres)
+{
+	uint32_t i, k, j, n, max, o;
+	uint8_t *mark = NULL;
+	CALLOC(mark, ma->ma.n);
+	for (i = 0; i < ma->n_seq; ++i) {
+		o = ma->idx.a[i] >> 32;
+		n = (uint32_t)ma->idx.a[i];
+		max = 0;
+		if (n == 0) continue;
+		for (j = o; j < o + n; ++j)
+			max = max > ma->ma.a[j].w? max : ma->ma.a[j].w;
+		for (j = o; j < o + n; ++j)
+			if (ma->ma.a[j].w >= (max*drop_thres) || ma->ma.a[j].w + min_cnt >= max)
+				mark[j] = 1;
+	}
+	for (i = 0; i < ma->ma.n; ++i)
+	{
+		if (mark[i] == 0) continue;
+		o = ma->idx.a[ma->ma.a[i].sid[1]]>>32;
+		n = (uint32_t)ma->idx.a[ma->ma.a[i].sid[1]];
+		for (j = o; j < o + n; ++j)
+		{
+			if (ma->ma.a[j].sid[1] == ma->ma.a[i].sid[0]) mark[j] = 1;
+		}
+	}
+		
+	for (i = k = 0; i < ma->ma.n; ++i)
+		if (mark[i]) ma->ma.a[k++] = ma->ma.a[i];
+	ma->ma.n = k;
+	free(mark);
+	pt_pdist_idx(ma);
+	pt_pdist_symm(ma);
+}
+
+pt_match_t *init_pt_match_t(hap_overlaps_list* ha, pt_g_t *x, double f_rate)
+{
+	pt_match_t *p = NULL; CALLOC(p, 1); p->n_seq = x->ug->g->n_seq;
+	kv_init(p->idx); kv_init(p->ma);
+	uint32_t v, i, k, qn, tn, qs, qe, ts, te, occ, as, ae;
+	uint64_t hetLen, homLen, oLen;
+	pt_node_t *a = NULL;
+	pt_match1_t *ma = NULL;
+
+	for (v = 0; v < ha->num; v++)
+    {
+        for (i = 0; i < ha->x[v].a.n; i++)
+        {
+			if(ha->x[v].a.a[i].score <= 0) continue;
+			if(ha->x[v].a.a[i].xUid == ha->x[v].a.a[i].yUid) continue;
+            /*****************qn*****************/
+            qn = ha->x[v].a.a[i].xUid;
+            qs = ha->x[v].a.a[i].x_beg_pos;
+            qe = ha->x[v].a.a[i].x_end_pos - 1;
+
+			a = x->p.a + x->p_idx.a[qn]; 
+			occ = x->p_idx.a[qn+1] - x->p_idx.a[qn];
+            for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
+            {
+                as = a[k].bS;
+                ae = a[k].bE;
+                oLen = ((MIN(qe, ae) >= MAX(qs, as))? MIN(qe, ae) - MAX(qs, as) + 1 : 0);
+                if(homLen + hetLen > 0 && oLen == 0) break;
+                if(oLen == 0) continue;
+                if(a[k].hs == N_HET)
+                {
+                    homLen += oLen;
+                } 
+                else if(asm_opt.polyploidy <= 2 && (a[k].hs&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].hs&S_HET))
+                {
+                    homLen += oLen;
+                }
+                else
+                {
+                    hetLen += oLen;
+                }                
+            }
+
+            if(hetLen <= ((hetLen + homLen)*f_rate)) continue;
+            /*****************qn*****************/
+
+
+            /*****************tn*****************/
+            tn = ha->x[v].a.a[i].yUid;
+            ts = ha->x[v].a.a[i].y_beg_pos;
+            te = ha->x[v].a.a[i].y_end_pos - 1;
+
+			a = x->p.a + x->p_idx.a[tn]; 
+			occ = x->p_idx.a[tn+1] - x->p_idx.a[tn];
+            for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
+            {
+                as = a[k].bS;
+                ae = a[k].bE;
+                oLen = ((MIN(te, ae) >= MAX(ts, as))? MIN(te, ae) - MAX(ts, as) + 1 : 0);
+                if(homLen + hetLen > 0 && oLen == 0) break;
+                if(oLen == 0) continue;
+                if(a[k].hs == N_HET)
+                {
+                    homLen += oLen;
+                } 
+                else if(asm_opt.polyploidy <= 2 && (a[k].hs&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].hs&S_HET))
+                {
+                    homLen += oLen;
+                }
+                else
+                {
+                    hetLen += oLen;
+                }                 
+            }
+
+            if(hetLen <= ((hetLen + homLen)*f_rate)) continue;
+            /*****************tn*****************/
+
+			kv_pushp(pt_match1_t, p->ma, &ma);
+			ma->sid[0] = ha->x[v].a.a[i].xUid;
+			ma->sid[1] = ha->x[v].a.a[i].yUid;
+			ma->w = ha->x[v].a.a[i].score;
+        }
+    }
+
+	pt_pdist_idx(p);
+	pt_pdist_symm(p);
+	pt_pdist_flt(p, 5, 0.5);
+	return p;
+}
+
+
+void debug_pt_g_t(pt_g_t *pg)
+{
+    fprintf(stderr, "0----------[M::%s]----------\n", __func__);
+    uint32_t i, offset, v, sid, eid, spos, epos, p_status, p_uid, occ;
+    ma_utg_t *u = NULL;
+    pt_node_t *a = NULL, *t = NULL;
+
+    for (v = 0; v < pg->ug->u.n; v++)
+    {
+		a = pg->p.a + pg->p_idx.a[v]; 
+		occ = pg->p_idx.a[v+1] - pg->p_idx.a[v];
+        for (i = 0; i < occ; i++)
+        {
+            if(a[i].uID != v) fprintf(stderr, "ERROR-s\n");
+        }
+    }
+
+    for (v = 0, p_status = (uint32_t)-1, p_uid = (uint32_t)-1; v < pg->p.n; v++)
+    {
+        t = &(pg->p.a[v]);
+        sid = t->nS;
+        eid = t->nE;
+        spos = t->bS;
+        epos = t->bE;
+        if(p_uid == t->uID && p_status == t->hs)
+        {
+            fprintf(stderr, "ERROR-a\n");
+        } 
+        p_status = t->hs;
+        p_uid = t->uID;
+        u = &(pg->ug->u.a[t->uID]);
+        for (i = offset = 0; i < u->n; i++)
+        {
+            if(i == sid)
+            {
+                if(spos != offset)
+                {
+                    fprintf(stderr, "ERROR-b\n");
+                }
+            }
+
+            if(i == eid)
+            {
+                if(epos != (offset+pg->rg->seq[u->a[i]>>33].len - 1))
+                {
+                    fprintf(stderr, "ERROR-c, real end: %u\n", 
+                                    (uint32_t)(offset+pg->rg->seq[u->a[i]>>33].len - 1));
+                }
+            }
+
+            offset += (uint32_t)u->a[i];
+            if(i >= sid && i <= eid)
+            {
+                if(pg->t_ch->is_r_het[u->a[i]>>33] != t->hs)
+                {
+                    fprintf(stderr, "ERROR-d: is_r_het: %u, h_status: %u\n", pg->t_ch->is_r_het[u->a[i]>>33], t->hs);
+                }
+            }
+        }
+    }
+
+
+	pt_match1_t *o = NULL, *s = NULL;
+	uint32_t k, n, cnt;
+	for (i = 0; i < pg->e->n_seq; ++i) 
+	{
+		o = pt_a(*(pg->e), i); n = pt_n(*(pg->e), i);
+		for (k = 0; k < n; ++k)
+		{
+			if(o[k].sid[0] != i) fprintf(stderr, "ERROR-g\n");
+			s = pt_pdist(pg->e, o[k].sid[1], o[k].sid[0]);
+			if(!s) fprintf(stderr, "ERROR-e\n");
+			if(s)
+			{
+				if(!(s->sid[0] == o[k].sid[1] && s->sid[1] == o[k].sid[0] && s->w == o[k].w))
+				{
+					fprintf(stderr, "ERROR-f\n");
+				}
+			}
+		}
+
+		for (k = cnt = 0; k < pg->e->ma.n; ++k)
+		{
+			if(pg->e->ma.a[k].sid[0] == i) cnt++;
+		}
+
+		if(cnt != n) fprintf(stderr, "ERROR-h\n");
+	}
+}
+
+pt_g_t *init_pt_g_t(hap_overlaps_list* ovlp, trans_chain* t_ch, ma_ug_t *ug, asg_t *read_g, double f_rate)
+{	
+	uint32_t v, l, k, offset, l_pos;
+	pt_g_t *p = NULL; CALLOC(p, 1);
+	asg_t* nsg = ug->g;
+	ma_utg_t *u = NULL;
+	pt_node_t *t = NULL;
+	p->ug = ug;
+	p->rg = read_g;
+	p->t_ch = t_ch;
+	kv_init(p->info); p->info.n = p->info.m = p->ug->g->n_seq; CALLOC(p->info.a, p->info.n);
+	kv_init(p->p);
+	kv_init(p->p_idx); kv_push(uint32_t, p->p_idx, 0);
+
+	for (v = 0; v < nsg->n_seq; v++)
+    {
+        if(nsg->seq[v].del || nsg->seq[v].c == ALTER_LABLE) continue;
+
+        u = &(ug->u.a[v]);
+        for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k) 
+        {   
+            if (k == u->n || t_ch->is_r_het[u->a[k]>>33] != t_ch->is_r_het[u->a[l]>>33])
+            {
+                kv_pushp(pt_node_t, p->p, &t);
+				t->uID = v;
+				t->hs = t_ch->is_r_het[u->a[l]>>33];
+
+				t->bS = l_pos;
+				t->bE = offset + read_g->seq[u->a[k-1]>>33].len - 1;
+
+				t->nS = l;
+				t->nE = k - 1;
+                l = k;
+                l_pos = offset + (uint32_t)u->a[k-1];
+            }
+            offset += (uint32_t)u->a[k-1];
+        }
+		kv_push(uint32_t, p->p_idx, p->p.n);
+    }
+
+	p->e = init_pt_match_t(ovlp, p, f_rate);
+	return p;
+}
+
+void destory_pt_g_t(pt_g_t **p)
+{
+	if(!p || !(*p)) return;
+	kv_destroy((*p)->p);
+	kv_destroy((*p)->info);
+	kv_destroy((*p)->p_idx);
+	kv_destroy((*p)->e->idx);
+	kv_destroy((*p)->e->ma);
+	free((*p)->e->cc);
+	free((*p)->e);
+	free((*p));
+}
+
+uint64_t *pt_cc_core(const pt_match_t *ma)
+{
+	uint32_t i, x, y, *flag;
+	uint64_t *group;
+	pt_match1_t *o = NULL;
+	kvec_t(uint32_t) stack; kv_init(stack);
+
+	MALLOC(flag, ma->n_seq);
+	for (i = 0; i < ma->n_seq; ++i)
+		flag[i] = (uint32_t)-1;
+
+	// connected componets
+	for (i = 0; i < ma->n_seq; ++i) {
+		if (flag[i] != (uint32_t)-1) continue;
+		stack.n = 0;
+		kv_push(uint32_t, stack, i);
+		while (stack.n > 0) {
+			uint32_t k, j, n;
+			stack.n--;
+			k = stack.a[stack.n];
+			flag[k] = i;///group id
+			// n = (uint32_t)ma->idx[k];
+			// s = ma->idx[k] >> 32;
+			o = pt_a(*ma, k); 
+			n = pt_n(*ma, k);
+			for (j = 0; j < n; ++j) {
+				uint32_t t = o[j].sid[1];
+				if (flag[t] != (uint32_t)-1) continue;
+				// if (ns == ms) PT_EXPAND(stack, ms);
+				// stack[ns++] = t;
+				kv_push(uint32_t, stack, t);
+			}
+		}
+	}
+	kv_destroy(stack);
+
+	// precalculate the size of each group
+	CALLOC(group, ma->n_seq);
+	for (i = 0; i < ma->n_seq; ++i)
+		group[i] = (uint64_t)flag[i] << 32 | i;
+	radix_sort_gfa64(group, group + ma->n_seq);
+	for (i = 1, x = y = 0; i <= ma->n_seq; ++i) {
+		if (i == ma->n_seq || group[i]>>32 != group[x]>>32) {
+			uint32_t j;
+			for (j = x; j < i; ++j)
+				group[j] = (uint64_t)y << 32 | (uint32_t)group[j];///(group id)|first element in this group
+			++y, x = i;
+		}
+	}
+	free(flag);
+	return group;
+}
+
+void pt_cc(pt_match_t *ma)
+{
+	ma->cc = pt_cc_core(ma);
+}
+
+///inspect top INT edges
+static int64_t pt_score(const pt_match_t *ma, uint32_t topn, solve_aux_t *aux)
+{
+	uint32_t i;
+	int64_t z = 0;
+	for (i = 0; i < aux->size; ++i) {///aux->size: how many unitigs in this group
+		uint32_t k = (uint32_t)ma->cc[aux->off + i];
+		uint32_t o = ma->idx.a[k] >> 32;
+		uint32_t n = (uint32_t)ma->idx.a[k], j;
+		for (j = 0; j < n; ++j)
+			aux->buf[j] = (uint64_t)((uint32_t)-1 - ma->ma.a[o + j].w) << 32 | (o + j);
+		radix_sort_gfa64(aux->buf, aux->buf + n);
+		for (j = 0; j < n && j < topn; ++j) {
+			const pt_match1_t *m = &ma->ma.a[(uint32_t)aux->buf[j]];
+			z += -(int64_t)m->w * aux->s[m->sid[0]] * aux->s[m->sid[1]];
+		}
+	}
+	return z;
+}
+
+static int64_t pt_solve1_init_phase(const pt_match_t *ma, int32_t topn, uint64_t *x, solve_aux_t *aux)
+{
+	uint32_t i;
+	aux->a.n = 0;
+	for (i = 0; i < aux->size; ++i) {
+		uint32_t k = (uint32_t)ma->cc[aux->off + i];///unitig id
+		uint32_t o = ma->idx.a[k] >> 32;///group id
+		uint32_t n = (uint32_t)ma->idx.a[k], j;
+		aux->shuffled[i] = k;///init shuffled status
+		for (j = 0; j < n; ++j) {
+			const pt_match1_t *m = &ma->ma.a[o + j];
+			///first is weight, second part is edge id
+			kv_push(uint64_t, aux->a, (uint64_t)((uint32_t)-1 - m->w) << 32 | (o + j));
+		}
+	}
+	radix_sort_gfa64(aux->a.a, aux->a.a + aux->a.n);///sort all edges in this group
+	///randomly assign each unitig
+	for (i = 0; i < aux->a.n; ++i) { // from the strongest edge to the weakest
+		const pt_match1_t *m = &ma->ma.a[(uint32_t)aux->a.a[i]];
+		///sid[0]: query id, sid[1]: target id
+		///the initial results of aux->s is 0
+		if (aux->s[m->sid[0]] == 0 && aux->s[m->sid[1]] == 0) {
+			*x = kr_splitmix64(*x);// random number
+			aux->s[m->sid[0]] = *x&1? 1 : -1;
+			aux->s[m->sid[1]] = -aux->s[m->sid[0]];
+		} else if (aux->s[m->sid[0]] == 0) {
+			aux->s[m->sid[0]] = -aux->s[m->sid[1]];
+		} else if (aux->s[m->sid[1]] == 0) {
+			aux->s[m->sid[1]] = -aux->s[m->sid[0]];
+		}
+	}
+	return pt_score(ma, topn, aux);
+}
+
+///size is how many unitigs in this group
+static void ks_shuffle_uint32_t(size_t n, uint32_t a[], uint64_t *x)
+{
+	size_t i, j;
+	for (i = n; i > 1; --i) {
+		uint32_t tmp;
+		j = (size_t)(kr_drand_r(x) * i);///semms 
+		tmp = a[j]; a[j] = a[i-1]; a[i-1] = tmp;
+	}
+}
+
+static void pt_solve1_perturb(const pt_svopt_t *opt, const pt_match_t *ma, uint64_t *x, solve_aux_t *aux)
+{
+	uint32_t i;
+	double y;
+	for (i = 0; i < aux->size; ++i) {
+		uint32_t k = (uint32_t)ma->cc[aux->off + i];
+		y = kr_drand_r(x);
+		if (y < opt->f_perturb)
+			aux->s[k] = -aux->s[k];
+	}
+}
+
+static int64_t pt_solve1_optimize(const pt_match_t *ma, uint32_t topn, uint64_t *x, solve_aux_t *aux, uint32_t *n_iter)
+{
+	uint32_t i;
+	while (1) {
+		uint32_t n_flip = 0;
+		++(*n_iter);
+		ks_shuffle_uint32_t(aux->size, aux->shuffled, x);
+		for (i = 0; i < aux->size; ++i) {
+			uint32_t k = aux->shuffled[i];
+			uint32_t o = ma->idx.a[k] >> 32;
+			uint32_t n = (uint32_t)ma->idx.a[k], j;
+			uint64_t z[2];
+			int8_t s;
+			for (j = 0; j < n; ++j) {
+				const pt_match1_t *m = &ma->ma.a[o + j];
+				///assert(m->sid[0] == k);
+				aux->buf[j] = (uint64_t)((uint32_t)-1 - m->w) << 32 | (o + j);
+			}
+			radix_sort_gfa64(aux->buf, aux->buf + n);///still sort by edge weight
+			for (j = 0, z[0] = z[1] = 0; j < n && j < topn; ++j) {
+				const pt_match1_t *m = &ma->ma.a[(uint32_t)aux->buf[j]];
+				if (aux->s[m->sid[1]] > 0) z[0] += m->w;
+				else if (aux->s[m->sid[1]] < 0) z[1] += m->w;
+			}
+			if (z[0] == z[1]) continue;
+			s = z[0] > z[1]? -1 : 1;
+			if (aux->s[k] != s)
+				aux->s[k] = s, ++n_flip;
+		}
+		if (n_flip == 0) break;
+	}
+	return pt_score(ma, topn, aux);
+}
+
+uint32_t pt_solve1(const pt_svopt_t *opt, const pt_match_t *ma, uint64_t *x, solve_aux_t *aux)
+{
+	uint32_t j, k, n_iter = 0;
+	int64_t sc_ori, sc_opt = -(1<<30), sc;
+	if (aux->size < 2) return 0;///how many unitigs
+
+	// first guess
+	///randomly assign haplotype status, and get a score
+	sc_ori = pt_solve1_init_phase(ma, opt->topn, x, aux);
+	if (aux->size == 2) return 0;
+
+	// optimize
+	sc_opt = pt_solve1_optimize(ma, opt->topn, x, aux, &n_iter);
+	for (j = 0; j < aux->size; ++j)
+		aux->s_tmp[aux->shuffled[j]] = aux->s[aux->shuffled[j]];
+	for (k = 0; k < (uint32_t)opt->n_perturb; ++k) {
+		pt_solve1_perturb(opt, ma, x, aux);
+		sc = pt_solve1_optimize(ma, opt->topn, x, aux, &n_iter);
+		if (sc > sc_opt) {
+			for (j = 0; j < aux->size; ++j)
+				aux->s_tmp[aux->shuffled[j]] = aux->s[aux->shuffled[j]];
+			sc_opt = sc;
+		} else {
+			for (j = 0; j < aux->size; ++j)
+				aux->s[aux->shuffled[j]] = aux->s_tmp[aux->shuffled[j]];
+		}
+	}
+	for (j = 0; j < aux->size; ++j)
+		aux->s[aux->shuffled[j]] = aux->s_tmp[aux->shuffled[j]];
+	fprintf(stderr, "[%s] group:%d, size:%d, #edges:%u, #iter:%d, sc_ori:%ld, sc_opt:%ld\n", __func__,
+			(uint32_t)(ma->cc[aux->off]>>32), (uint32_t)(aux->size), (uint32_t)(aux->a.n), n_iter, (long)sc_ori, (long)sc_opt);
+	return n_iter;
+}
+
+int8_t *pt_solve_core(const pt_svopt_t *opt, const pt_match_t *ma)
+{
+	int8_t *s;
+	uint32_t st, i, max = 0;
+	uint64_t x = opt->seed;
+	solve_aux_t *aux;
+	CALLOC(aux, 1); CALLOC(aux->s, ma->n_seq); CALLOC(aux->s_tmp, ma->n_seq);
+	kv_init(aux->a);
+	for (i = 0; i < ma->n_seq; ++i) {///count how many links for each unitig
+		uint32_t n = pt_n(*ma, i);
+		max = max > n? max : n;
+	}
+	MALLOC(aux->buf, max);
+	MALLOC(aux->shuffled, ma->n_seq); // FIXME: this is over-allocation for convenience
+	for (st = 0, i = 1; i <= ma->n_seq; ++i) {
+		if (i == ma->n_seq || ma->cc[st]>>32 != ma->cc[i]>>32) {///at same group
+			if (i - st >= 2) {///all unitigs in the same group
+				aux->off = st, aux->size = i - st;
+				pt_solve1(opt, ma, &x, aux);
+			}
+			st = i;
+		}
+	}
+	s = aux->s;
+	kv_destroy(aux->a); free(aux->buf); free(aux->shuffled); free(aux->s_tmp);
+	free(aux);
+	return s;
+}
+
+void set_trio_flag(ma_ug_t *ug, asg_t *read_g, uint32_t uID, uint8_t* trio_flag, trans_chain* t_ch, 
+													hap_overlaps_list* ha, pt_match_t *ma, int8_t s)
+{
+	uint32_t i;
+	ma_utg_t *u = &(ug->u.a[uID]);
+	for (i = 0; i < u->n; i++)
+	{
+		trio_flag[u->a[i]>>33] |= SET_TRIO;
+		if(t_ch->is_r_het[u->a[i]>>33] == N_HET) continue;
+		if(s == 0)
+		{
+			if(t_ch->is_r_het[u->a[i]>>33]&P_HET)//special case
+			{
+				trio_flag[u->a[i]>>33] |= FATHER;
+			}
+			continue;
+		} 
+		trio_flag[u->a[i]>>33] |= (s > 0? FATHER:MOTHER);
+	}
+
+	// pt_match1_t *o = pt_a(*ma, uID);
+	// uint32_t n = pt_n(*ma, uID), k, qn, tn, qs, qe, r_beg, r_end, offset, oLen, found;
+	// int index;
+	// for (k = 0; k < n; ++k)
+	// {
+	// 	qn = o[k].sid[0]; tn = o[k].sid[1]; qs = qe = (uint32_t)-1;
+	// 	index = get_specific_hap_overlap(&(ha->x[qn]), qn, tn);
+	// 	if(index != -1 && ha->x[qn].a.a[index].score == (long long)o[k].w)
+	// 	{
+	// 		qs = ha->x[qn].a.a[index].x_beg_pos;
+    //         qe = ha->x[qn].a.a[index].x_end_pos - 1;
+	// 	}
+	// 	else
+	// 	{
+	// 		index = get_specific_hap_overlap(&(ha->x[tn]), tn, qn);
+	// 		if(index != -1 && ha->x[tn].a.a[index].score == (long long)o[k].w)
+	// 		{
+	// 			qs = ha->x[qn].a.a[index].y_beg_pos;
+    //         	qe = ha->x[qn].a.a[index].y_end_pos - 1;
+	// 		}
+	// 	}
+		
+	// 	if(qs == (uint32_t)-1 || qe == (uint32_t)-1) fprintf(stderr, "ERROR\n");
+	// 	for (i = 0, offset = 0, found = 0; i < u->n; i++)
+	// 	{
+	// 		r_beg = offset; r_end = offset + (long long)(read_g->seq[u->a[i]>>33].len) - 1;
+    //     	offset += (uint32_t)u->a[i];
+	// 		oLen = ((MIN(qe, r_end) >= MAX(qs, r_beg))? MIN(qe, r_end) - MAX(qs, r_beg) + 1 : 0);
+	// 		if(oLen > 0) found = 1;
+	// 		if(found == 1 && oLen == 0) break;
+
+	// 		if(oLen > 0 && t_ch->is_r_het[u->a[i]>>33] != N_HET)
+	// 		{
+	// 			trio_flag[u->a[i]>>33] |= (s > 0? FATHER:MOTHER);
+	// 		} 
+			
+	// 	}
+	// }
+}
+
+void filter_ovlp(ma_ug_t *ug, asg_t *read_g, uint32_t uID, hap_overlaps_list* ha, pt_match_t *ma, 
+int8_t *s)
+{
+	pt_match1_t *o = pt_a(*ma, uID);
+	uint32_t n = pt_n(*ma, uID), k, qn, tn;
+	hap_overlaps *p = NULL;
+	int index;
+
+	for (k = 0; k < n; ++k)
+	{
+		qn = o[k].sid[0]; tn = o[k].sid[1]; p = NULL;
+		if((s[qn]*s[tn])!=-1) continue;
+		index = get_specific_hap_overlap(&(ha->x[qn]), qn, tn);
+		if(index != -1 && ha->x[qn].a.a[index].score == (long long)o[k].w)
+		{
+			p = &(ha->x[qn].a.a[index]);
+		}
+		else
+		{
+			index = get_specific_hap_overlap(&(ha->x[tn]), tn, qn);
+			if(index != -1 && ha->x[tn].a.a[index].score == (long long)o[k].w)
+			{
+				p = &(ha->x[tn].a.a[index]);
+			}
+		}
+		if(!p) fprintf(stderr, "ERROR\n");
+		p->status = FLIP;
+	}
+}
+
+
+void clean_ovlp(ma_ug_t *ug, asg_t *read_g, hap_overlaps_list* ha, pt_g_t *pg, int8_t* s)
+{
+	uint32_t v, i, k, qn, tn, types[4];
+    types[X2Y] = Y2X; types[Y2X] = X2Y; types[XCY] = YCX; types[YCX] = XCY;
+	int index;
+	hap_overlaps *x = NULL, *y = NULL;
+	for (i = 0; i < pg->e->n_seq; ++i) 
+	{
+		filter_ovlp(ug, read_g, i, ha, pg->e, s);
+	}
+
+	for (v = 0; v < ha->num; v++)
+	{
+		for (i = 0; i < ha->x[v].a.n; i++)
+		{
+			qn = ha->x[v].a.a[i].xUid;
+			tn = ha->x[v].a.a[i].yUid;
+			x = &(ha->x[v].a.a[i]);
+			if(x->status != FLIP) continue;
+			index = get_specific_hap_overlap(&(ha->x[tn]), tn, qn);
+			if(index != -1)
+			{
+				y = &(ha->x[tn].a.a[index]);
+				set_reverse_hap_overlap(y, x, types);
+				y->status = FLIP;
+			}
+			if(index == -1) fprintf(stderr, "ERROR\n");
+		}
+	}
+
+
+	// for (v = 0; v < ha->num; v++)
+    // {
+    //     for (i = k = 0; i < ha->x[v].a.n; i++)
+    //     {
+    //         if(ha->x[v].a.a[i].status != FLIP)
+	// 		{
+	// 			if(s[ha->x[v].a.a[i].xUid]*s[ha->x[v].a.a[i].yUid] == -1)
+	// 			{
+	// 				fprintf(stderr, "\ns[0]=%d, s[1]=%d\n", s[ha->x[v].a.a[i].xUid], s[ha->x[v].a.a[i].yUid]);
+	// 				print_hap_paf(ug, &(ha->x[v].a.a[i]));
+	// 			}
+				
+	// 			continue;
+	// 		} 
+			
+			
+    //         ha->x[v].a.a[k] = ha->x[v].a.a[i];
+    //         k++;
+    //     }
+    //     ha->x[v].a.n = k;
+    // }
+
+	for (v = 0; v < ha->num; v++)
+    {
+        for (i = k = 0; i < ha->x[v].a.n; i++)
+        {
+            if(ha->x[v].a.a[i].status != FLIP) continue;
+
+            ha->x[v].a.a[k] = ha->x[v].a.a[i];
+            k++;
+        }
+        ha->x[v].a.n = k;
+    }
+}
+
+void pt_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag)
+{
+	pt_svopt_t opt;
+	int8_t *s = NULL;
+	uint64_t *buf = NULL, i;
+	pt_svopt_init(&opt);
+	pt_g_t *pg = init_pt_g_t(ovlp, t_ch, ug, read_g, f_rate);
+	///debug_pt_g_t(pg);
+
+	pt_cc(pg->e);
+	s = pt_solve_core(&opt, pg->e);
+
+	if(asm_opt.flag & HA_F_PARTITION)
+	{
+		MALLOC(buf, pg->e->ma.n); // FIXME: this is over-allocation for convenience
+		for (i = 0; i < pg->e->n_seq; ++i) {
+			uint64_t z[2];
+			uint32_t o = pg->e->idx.a[i] >> 32;
+			uint32_t n = (uint32_t)pg->e->idx.a[i], j;
+
+			set_trio_flag(ug, read_g, i, trio_flag, t_ch, ovlp, pg->e, s[i]);
+
+			pg->info.a[i].s = s[i];
+			for (j = 0; j < n; ++j) {
+				const pt_match1_t *m = &pg->e->ma.a[o + j];
+				buf[j] = (uint64_t)((uint32_t)-1 - m->w) << 32 | (o + j);
+			}
+			radix_sort_gfa64(buf, buf + n);
+			for (j = 0, z[0] = z[1] = 0; j < n; ++j) {
+				const pt_match1_t *m = &pg->e->ma.a[(uint32_t)buf[j]];
+				if (s[m->sid[1]] > 0) z[0] += m->w;
+				else if (s[m->sid[1]] < 0) z[1] += m->w;
+			}
+			pg->info.a[i].m[0] = z[0], pg->info.a[i].m[1] = z[1];
+		}
+	}
+	
+
+	clean_ovlp(ug, read_g, ovlp, pg, s);
+
+	free(buf);
+	free(s);
+	destory_pt_g_t(&pg);
+}

partig.h

@@ -0,0 +1,51 @@
+#ifndef __PURTIG__
+#define __PURTIG__
+#include <stdio.h>
+#include <stdint.h>
+#include "kvec.h"
+#include "Overlaps.h"
+#include "Purge_Dups.h"
+
+typedef struct {
+    uint32_t bS, bE;
+    uint32_t nS, nE;
+    uint32_t uID;
+    uint8_t hs;
+}pt_node_t;
+
+typedef struct {
+	///sid[0]: query id
+	///sid[1]: target id
+	uint32_t sid[2];
+	uint32_t w;
+} pt_match1_t;
+
+typedef struct {
+    kvec_t(uint64_t) idx;
+    kvec_t(pt_match1_t) ma;
+    uint64_t* cc;
+    uint32_t n_seq;
+} pt_match_t;
+
+typedef struct {
+	///cnt1: how many unique minimizers
+	///cnt2: how many non-unique minimizers
+	///uint32_t cnt2, cnt1;
+	uint64_t m[2];
+	int8_t s;
+} pt_uinfo_t;
+
+typedef struct {
+    kvec_t(pt_node_t) p;
+    ma_ug_t *ug;
+    asg_t *rg;
+    trans_chain* t_ch;
+    ///kvec_t(int8_t) s; ///status
+    kvec_t(pt_uinfo_t) info; ///status
+    kvec_t(uint32_t) p_idx;
+    pt_match_t* e;
+}pt_g_t;
+
+void pt_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag);
+
+#endif

rcut.h

@@ -0,0 +1,62 @@
+#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 int32_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)))
+
+
+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;
+
+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);
+#endif