fix version number

Hifiasm high het

CommandLines.cpp

@@ -78,7 +78,7 @@ void Print_H(hifiasm_opt_t* asm_opt)
                                      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\n");
+    fprintf(stderr, "    --high-het  enable this mode for high heterozygosity sample [experimental, not stable]\n");
     
 
     fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.9-r297"
+#define HA_VERSION "0.11-r302"
 
 #define VERBOSE 0
 

Makefile

@@ -1,5 +1,7 @@
 CXX=		g++
+CC=			gcc
 CXXFLAGS=	-g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall
+CFLAGS=		$(CXXFLAGS)
 CPPFLAGS=
 INCLUDES=
 OBJS=		CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
@@ -13,12 +15,15 @@ ifneq ($(asan),)
 	LIBS+=-fsanitize=address
 endif
 
-.SUFFIXES:.cpp .o
+.SUFFIXES:.cpp .c .o
 .PHONY:all clean depend
 
 .cpp.o:
 		$(CXX) -c $(CXXFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@
 
+.c.o:
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@
+
 all:$(EXE)
 
 $(EXE):$(OBJS) main.o
@@ -36,8 +41,9 @@ Assembly.o: Assembly.h CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h
 Assembly.o: Hash_Table.h htab.h POA.h Correct.h Levenshtein_distance.h
 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 ksw2.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
 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
@@ -65,4 +71,3 @@ main.o: CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h Assembly.h
 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
-ksw2_extz2_sse.o: ksw2.h

Overlaps.cpp

@@ -6618,7 +6618,7 @@ uint8_t get_tip_trio_infor(asg_t *sg, uint32_t begNode)
 {
     uint32_t v = begNode, w;
     uint32_t kv;
-    uint32_t eLen = 0;
+    uint32_t eLen = 0, uLen = 0;
     uint32_t father_occ = 0, mother_occ = 0, ambigious_occ = 0;
 
     while (1)
@@ -6647,10 +6647,17 @@ uint8_t get_tip_trio_infor(asg_t *sg, uint32_t begNode)
         if(v == begNode) break;
     }
 
+    uLen = eLen;
     eLen = father_occ + mother_occ;
-    if(father_occ >= TRIO_THRES*eLen) return FATHER;
+    if(eLen == 0) return AMBIGU;
-    if(mother_occ >= TRIO_THRES*eLen) return MOTHER;
+    if(father_occ >= mother_occ)
-    
+    {
+        if((father_occ > TRIO_THRES*eLen) && (father_occ >= DOUBLE_CHECK_THRES*uLen)) return FATHER;
+    }
+    else
+    {
+        if((mother_occ > TRIO_THRES*eLen) && (mother_occ >= DOUBLE_CHECK_THRES*uLen)) return MOTHER;
+    }
     return AMBIGU;
 }
 
@@ -6697,28 +6704,30 @@ int if_skip_bubble, int if_drop, int if_check_hap, R_to_U* ruIndex)
         kv = get_real_length(g, v, NULL);
         kw = get_real_length(g, w, NULL);
         if (kv <= 1 && kw <= 1) continue;
-        trio_flag = get_tip_trio_infor(g, v);
+        trio_flag = get_tip_trio_infor(g, v^1);
         non_trio_flag = (uint32_t)-1;
         if(trio_flag == FATHER) non_trio_flag = MOTHER;
         if(trio_flag == MOTHER) non_trio_flag = FATHER;
-
+        
         ///calculate the longest edge for v and w
 		for (i = 0, kv = 0; i < nv; ++i) {
 			if (av[i].del) continue;
+            kv++; 
             if(get_tip_trio_infor(g, av[i].v) == non_trio_flag) continue;
 			if (ov_max < av[i].ol) ov_max = av[i].ol, ov_max_i = i;
-            kv++; 
+            ///kv++; 
 		}
 		if (kv >= 2 && a->ol > ov_max * drop_ratio) continue;
 
-
 		for (i = 0, kw = 0; i < nw; ++i) {
 			if (aw[i].del) continue;
+            kw++;
             if (get_tip_trio_infor(g, aw[i].v) == non_trio_flag) continue;
 			if (ow_max < aw[i].ol) ow_max = aw[i].ol, ow_max_i = i;
-            kw++;
+            ///kw++;
 		}
 		if (kw >= 2 && a->ol > ow_max * drop_ratio) continue;
+
         if (kv <= 1 && kw <= 1) continue;
 
         ///to see which one is the current edge (from v and w)
@@ -7999,7 +8008,7 @@ int asg_arc_del_short_diploid_by_exact_trio(asg_t *g, int max_ext, ma_hit_t_allo
         kv = get_real_length(g, v, NULL);
         kw = get_real_length(g, w, NULL);
         if (kv <= 1 && kw <= 1) continue;
-        trio_flag = get_tip_trio_infor(g, v);
+        trio_flag = get_tip_trio_infor(g, v^1);
         non_trio_flag = (uint32_t)-1;
         if(trio_flag == FATHER) non_trio_flag = MOTHER;
         if(trio_flag == MOTHER) non_trio_flag = FATHER;
@@ -8008,25 +8017,27 @@ int asg_arc_del_short_diploid_by_exact_trio(asg_t *g, int max_ext, ma_hit_t_allo
         ///calculate the longest edge for v and w
 		for (i = 0, kv = 0; i < nv; ++i) {
 			if (av[i].del) continue;
+            ++kv;
             if(get_tip_trio_infor(g, av[i].v) == non_trio_flag) continue;
 			if (ov_max < av[i].ol) 
             {
                 ov_max = av[i].ol;
                 ov_max_i = i;
             }
-			++kv;
+			///++kv;
 		}
 		if (kv >= 2 && a->ol == ov_max) continue;
 
 
 		for (i = 0, kw = 0; i < nw; ++i) {
 			if (aw[i].del) continue;
+            ++kw;
             if (get_tip_trio_infor(g, aw[i].v) == non_trio_flag) continue;
 			if (ow_max < aw[i].ol) 
             {
                 ow_max = aw[i].ol;
             }
-			++kw;
+			///++kw;
 		}
 		if (kw >= 2 && a->ol == ow_max) continue;
         if (kv <= 1 && kw <= 1) continue;
@@ -9571,20 +9582,20 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
 
 
 
-void debug_info_of_specfic_node(char* name, asg_t *g, char* command)
+void debug_info_of_specfic_node(const char* name, asg_t *g, char* command)
 {
     fprintf(stderr, "\n\n\n");
-    uint32_t v, n_vtx = g->n_seq * 2;
+    uint32_t v, n_vtx = g->n_seq * 2, queryLen = strlen(name);
     for (v = 0; v < n_vtx; ++v) 
     {
-        if(memcmp(name, Get_NAME(R_INF, (v>>1)), Get_NAME_LENGTH(R_INF, (v>>1))) == 0)
+        if(queryLen == Get_NAME_LENGTH(R_INF, (v>>1)) && memcmp(name, Get_NAME(R_INF, (v>>1)), Get_NAME_LENGTH(R_INF, (v>>1))) == 0)
         {
             fprintf(stderr, "\nafter %s\n****************graph ref_read: %.*s, dir: %u****************\n", 
             command, (int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), v&1);
             if(g->seq[v>>1].del)
             {
                 fprintf(stderr, "read has already been deleted.\n");
-                continue;
+                return;
             }
 
             asg_arc_t *av = asg_arc_a(g, v);
@@ -9596,6 +9607,7 @@ void debug_info_of_specfic_node(char* name, asg_t *g, char* command)
                 Get_NAME(R_INF, (av[i].v>>1)),
                 av[i].el, av[i].strong, av[i].ol, av[i].del);
             }
+            return;
         }
     }
 }
@@ -11356,6 +11368,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp)
     kvec_asg_arc_t_warp new_rtg_edges;
     kv_init(new_rtg_edges.a);
 
+    if(ug == NULL) ug = ma_ug_gen(read_g);
 
     uint32_t i;
     for (i = 0; i < ug->u.n; ++i) 
@@ -26828,8 +26841,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     ///just need to deal with trio here
     ma_hit_contained_advance(sources, n_read, coverage_cut, ruIndex, max_hang_length, mini_overlap_length);
 
-    // debug_info_of_specfic_read("m54329U_190617_231905/176226420/ccs", sources, reverse_sources, 
+    ///debug_info_of_specfic_read("m54329U_190827_173812/166332272/ccs", sources, reverse_sources, -1, "clean");
-    // -1, "clean");
 
     sg = ma_sg_gen(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
     asg_arc_del_trans(sg, gap_fuzz);
@@ -26846,7 +26858,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     asg_cut_tip(sg, asm_opt.max_short_tip);
     
     ///drop_inexact_edegs_at_bubbles(sg, bubble_dist);
-
+    ///debug_info_of_specfic_node("m54329U_190827_173812/166332272/ccs", sg, "beg");
+    
     if(clean_round > 0)
     {
         double cut_step;
@@ -26862,7 +26875,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
         int i = 0;
         for (i = 0; i < clean_round; i++, drop_ratio += cut_step)
         {
-
             if(drop_ratio > max_ovlp_drop_ratio)
             {
                 drop_ratio = max_ovlp_drop_ratio;
@@ -26873,22 +26885,22 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
                 fprintf(stderr, "\n\n**********%d-th round drop: drop_ratio = %f**********\n", 
                 i, drop_ratio);
             }
-            
+
             ///just topological clean
             pre_clean(sources, coverage_cut, sg, bubble_dist);
-
             ///asg_arc_del_orthology(sg, reverse_sources, drop_ratio, asm_opt.max_short_tip);
             // asg_arc_del_orthology_multiple_way(sg, reverse_sources, drop_ratio, asm_opt.max_short_tip);
             // asg_cut_tip(sg, asm_opt.max_short_tip);
             /****************************may have bugs********************************/
+
             asg_arc_identify_simple_bubbles_multi(sg, 1);
             //reomve edge between two chromesomes
             //this node must be a single read
             asg_arc_del_false_node(sg, asm_opt.max_short_tip);
             asg_cut_tip(sg, asm_opt.max_short_tip);
             /****************************may have bugs********************************/
-
             /****************************may have bugs********************************/
+
             ///asg_arc_identify_simple_bubbles_multi(sg, 1);
             asg_arc_identify_simple_bubbles_multi(sg, 0);
             ///asg_arc_del_short_diploid_unclean_exact(sg, drop_ratio, sources);
@@ -26903,7 +26915,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
             asg_cut_tip(sg, asm_opt.max_short_tip);
             /****************************may have bugs********************************/
 
-
             asg_arc_identify_simple_bubbles_multi(sg, 1);
             if (ha_opt_triobin(&asm_opt))
             {
@@ -26919,7 +26930,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
 
             asg_arc_identify_simple_bubbles_multi(sg, 1);
             asg_arc_del_short_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, 
-            asm_opt.max_short_tip, ruIndex);
+            asm_opt.max_short_tip, ruIndex);     
 
             asg_arc_identify_simple_bubbles_multi(sg, 1);
             asg_arc_del_complex_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, asm_opt.max_short_tip);

hifiasm.1

@@ -264,6 +264,12 @@ Min number of overlapped reads for duplicate haplotigs that should be purged [1]
 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].
+
+
 .SS Debugging options
 
 .TP 10