Merge pull request #8 from chhylp123/flt-p

for purge_dup

CommandLines.h

@@ -4,7 +4,7 @@
 #include <pthread.h>
 
 #define VERBOSE 0
-#define VERBOSE_GFA 0
+#define VERBOSE_GFA 1
 
 typedef struct {
     int num_reads;

Hash_Table.cpp

@@ -672,7 +672,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
 
 
     clear_fake_cigar(&(result->f_cigar));
-    ///not a has been sorted by offset, that means has been sorted by query offset
+    ///note a has been sorted by offset, that means has been sorted by query offset
     i = max_i;
     result->x_pos_e = a[i].self_offset;
     result->y_pos_e = a[i].offset;

Hash_Table.h

@@ -549,6 +549,7 @@ void init_window_list_alloc(window_list_alloc* x);
 void clear_window_list_alloc(window_list_alloc* x);
 void destory_window_list_alloc(window_list_alloc* x);
 void resize_window_list_alloc(window_list_alloc* x, long long size);
-
+void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, 
+double band_width_threshold, int max_skip, int x_readLen, int y_readLen);
 
 #endif

Makefile

@@ -3,7 +3,7 @@ CXXFLAGS=	-g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall #-fsanitize=addres
 CPPFLAGS=
 INCLUDES=
 OBJS=		Output.o CommandLines.o Process_Read.o Assembly.o kmer.o Hash_Table.o \
-			POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o #ksw2_extz2_sse.o
+			POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o #ksw2_extz2_sse.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm #-fsanitize=address -fno-omit-frame-pointer
 
@@ -44,7 +44,7 @@ Hash_Table.o: POA.h ksort.h
 Levenshtein_distance.o: Levenshtein_distance.h
 Output.o: Output.h CommandLines.h
 Overlaps.o: Overlaps.h kvec.h kdq.h ksort.h Process_Read.h kseq.h
-Overlaps.o: CommandLines.h
+Overlaps.o: CommandLines.h Purge_Dups.h
 POA.o: POA.h Hash_Table.h khash.h kmer.h Process_Read.h kseq.h Overlaps.h
 POA.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h
 Process_Read.o: Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h
@@ -53,4 +53,5 @@ main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h
 main.o: Assembly.h Levenshtein_distance.h
 Trio.o: Trio.h khashl.h kthread.h Process_Read.h CommandLines.h
 kthread.o: kthread.h
+Purge_Dups.o: Purge_Dups.h Overlaps.h Hash_Table.h Correct.h ksort.h
 #ksw2_extz2_sse.o: ksw2.h

Overlaps.h

@@ -137,6 +137,16 @@ typedef struct {
 	ma_utg_t* F_seq;
 } asg_t;
 
+asg_t *asg_init(void);
+void asg_destroy(asg_t *g);
+void asg_arc_sort(asg_t *g);
+void asg_seq_set(asg_t *g, int sid, int len, int del);
+void asg_arc_index(asg_t *g);
+void asg_cleanup(asg_t *g);
+void asg_symm(asg_t *g);
+void print_gfa(asg_t *g);
+
+
 typedef struct { size_t n, m; uint64_t *a; } asg64_v;
 
 
@@ -373,6 +383,7 @@ typedef struct {
 	uint64_t i;
 }kvec_asg_arc_t_warp;
 
+void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend);
 
 
 typedef struct {
@@ -777,7 +788,10 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
         	{
 				qn = node_min->a[b_min.readI]>>33;
 
-				if(reverse_sources[qn].length>=0) min_count++;
+				/************************BUG: don't forget****************************/
+				if(reverse_sources[qn].length > 0) min_count++;
+				///if(reverse_sources[qn].length >= 0) min_count++;
+				/************************BUG: don't forget****************************/
 				for (j = 0; j < (long long)reverse_sources[qn].length; j++)
 				{
 					tn = Get_tn(reverse_sources[qn].buffer[j]);
@@ -826,7 +840,10 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 		{
 			qn = (b_min.b_0->b.a[b_min.untigI]>>1);
 
-			if(reverse_sources[qn].length>=0) min_count++;
+			/************************BUG: don't forget****************************/
+			if(reverse_sources[qn].length > 0) min_count++;
+			///if(reverse_sources[qn].length >= 0) min_count++;
+			/************************BUG: don't forget****************************/
 
 			for (j = 0; j < (long long)reverse_sources[qn].length; j++)
         	{
@@ -870,6 +887,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 }
 
 
+
 inline uint32_t check_different_haps_naive(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg, 
 uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1, 
 R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
@@ -923,7 +941,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
         	{
 				qn = node_min->a[b_min.readI]>>33;
 
-				if(reverse_sources[qn].length>=0) min_count++;
+				if(reverse_sources[qn].length > 0) min_count++;
 				for (j = 0; j < (long long)reverse_sources[qn].length; j++)
 				{
 					tn = Get_tn(reverse_sources[qn].buffer[j]);
@@ -960,7 +978,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 		{
 			qn = (b_min.b_0->b.a[b_min.untigI]>>1);
 
-			if(reverse_sources[qn].length>=0) min_count++;
+			if(reverse_sources[qn].length > 0) min_count++;
 
 			for (j = 0; j < (long long)reverse_sources[qn].length; j++)
         	{
@@ -1020,7 +1038,7 @@ void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_allo
 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);
 void deduplicate(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);
+long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t resolve_tangle);
 void all_to_all_deduplicate(ma_ug_t* ug, uint8_t postive_flag, float drop_rate, 
 ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
 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);
@@ -1032,5 +1050,11 @@ R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t bac
 uint32_t is_bubble_check, uint32_t is_primary_check);
 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);
+int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio);
+
+#define JUNK_COV 5
+#define DISCARD_RATE 0.8
 
 #endif

Purge_Dups.h

@@ -0,0 +1,19 @@
+#ifndef __PURGEDUPS__
+#define __PURGEDUPS__
+#include <stdio.h>
+#include <stdint.h>
+#include "kvec.h"
+#include "kdq.h"
+#include "Overlaps.h"
+#include "Hash_Table.h"
+
+void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* reverse_sources, 
+R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t bi_graph_Len, uint32_t long_hap_overlap, 
+float lable_match_rate, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio,
+uint32_t just_contain);
+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 enable_debug_mode();
+
+#endif