update for 0.2.0

Assembly.cpp

@@ -1716,7 +1716,7 @@ void Output_PAF()
 {
 
     fprintf(stderr, "Writing PAF to disk ...... \n");
-    char* paf_name = (char*)malloc(strlen(asm_opt.output_file_name)+5);
+    char* paf_name = (char*)malloc(strlen(asm_opt.output_file_name)+50);
     sprintf(paf_name, "%s.ovlp.paf", asm_opt.output_file_name);
     FILE* output_file = fopen(paf_name, "w");
     uint64_t i, j;
@@ -1757,6 +1757,8 @@ void Output_PAF()
 
     free(paf_name);
     fclose(output_file);
+
+    fprintf(stderr, "PAF has beem written.\n");
 }
 
 

CommandLines.cpp

@@ -5,7 +5,7 @@
 #include "ketopt.h"
 #include <sys/time.h>
 
-#define VERSION "0.1.0"
+#define VERSION "0.2.0"
 #define DEFAULT_OUTPUT "hifiasm.asm"
 
 hifiasm_opt_t asm_opt;

Correct.cpp

@@ -5024,7 +5024,8 @@ long long xBeg, long long xEnd, long long flag_offset)
             for (i = 0; i < operationLen; i++)
             {
                 /// note we need to deal with flag_offset carefully
-                if(flag[x_i - flag_offset] < 127 && x_i >= xBeg && x_i <= xEnd)
+                ///if(flag[x_i - flag_offset] < 127 && x_i >= xBeg && x_i <= xEnd)
+                if(x_i >= xBeg && x_i <= xEnd && flag[x_i - flag_offset] < 127)
                 {
                     flag[x_i - flag_offset]++;
                 }
@@ -5275,7 +5276,8 @@ haplotype_evdience_alloc* hap, long long snp_threshold)
             {
                 ///should be at least 2 mismatches
                 /// note we need to deal with flag_offset carefully
-                if(flag[x_i - flag_offset] > snp_threshold && x_i >= xBeg && x_i <= xEnd)
+                ///if(flag[x_i - flag_offset] > snp_threshold && x_i >= xBeg && x_i <= xEnd)
+                if(x_i >= xBeg && x_i <= xEnd && flag[x_i - flag_offset] > snp_threshold)
                 {
                     ev.misBase =  y_string[y_i];
                     ev.overlapID = overlapID;
@@ -5296,7 +5298,8 @@ haplotype_evdience_alloc* hap, long long snp_threshold)
 
                 /// should be at least 2 mismatches
                 /// note we need to deal with flag_offset carefully
-                if(flag[x_i - flag_offset] > snp_threshold && x_i >= xBeg && x_i <= xEnd)
+                ///if(flag[x_i - flag_offset] > snp_threshold && x_i >= xBeg && x_i <= xEnd)
+                if(x_i >= xBeg && x_i <= xEnd && flag[x_i - flag_offset] > snp_threshold)
                 {
                     ev.misBase =  y_string[y_i];
                     ev.overlapID = overlapID;
@@ -5323,7 +5326,8 @@ haplotype_evdience_alloc* hap, long long snp_threshold)
                 ///if(hap->flag[inner_offset] > snp_threshold)
                 /// should be at least 2 mismatches
                 /// note we need to deal with flag_offset carefully
-                if(flag[x_i - flag_offset] > snp_threshold && x_i >= xBeg && x_i <= xEnd)
+                ///if(flag[x_i - flag_offset] > snp_threshold && x_i >= xBeg && x_i <= xEnd)
+                if(x_i >= xBeg && x_i <= xEnd && flag[x_i - flag_offset] > snp_threshold)
                 {
                     ev.misBase =  'N';
                     ev.overlapID = overlapID;

Makefile

@@ -1,11 +1,11 @@
 CXX=		g++
-CXXFLAGS=	-g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall #-Winline 
+CXXFLAGS=	-g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall #-fsanitize=address -fno-omit-frame-pointer#-Winline 
 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 #ksw2_extz2_sse.o
 EXE=		hifiasm
-LIBS=		-lz -lpthread -lm
+LIBS=		-lz -lpthread -lm #-fsanitize=address -fno-omit-frame-pointer
 
 ifneq ($(asan),)
 	CXXFLAGS+=-fsanitize=address

Overlaps.h

@@ -17,7 +17,12 @@
 ///#define MAX_BUBBLE_DIST 10000000
 #define SMALL_BUBBLE_SIZE (uint32_t)-1
 //#define SMALL_BUBBLE_SIZE 1000
-
+#define PRIMARY_LABLE 0
+#define ALTER_LABLE 1
+#define HAP_LABLE 2
+// #define PRIMARY_LABLE 1
+// #define ALTER_LABLE 2
+// #define HAP_LABLE 4
 
 #define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32))
 #define Get_qs(RECORD) ((uint32_t)((RECORD).qns))
@@ -36,6 +41,8 @@
 #define LOOP 7
 
 
+
+
 ///query is the read itself
 typedef struct {
 	uint64_t qns;
@@ -46,7 +53,6 @@ typedef struct {
 	uint8_t no_l_indel;
 } ma_hit_t;
 
-
 typedef struct {
 	ma_hit_t* buffer;
     uint32_t size;
@@ -241,7 +247,7 @@ static inline void asg_seq_del(asg_t *g, uint32_t s)
 static inline void asg_seq_drop(asg_t *g, uint32_t s)
 {
 	///s is not at primary
-	if(g->seq[s].c)
+	if(g->seq[s].c == ALTER_LABLE)
 	{
 		uint32_t k;
 		for (k = 0; k < 2; ++k) 
@@ -254,8 +260,10 @@ static inline void asg_seq_drop(asg_t *g, uint32_t s)
 			{
 				if(av[i].del) continue;
 				///if output node is at primary
-				if(g->seq[(av[i].v>>1)].c == 0)
+				/****************************may have hap bugs********************************/
-				{
+				///if(g->seq[(av[i].v>>1)].c == PRIMARY_LABLE)
+				if(g->seq[(av[i].v>>1)].c == PRIMARY_LABLE || g->seq[(av[i].v>>1)].c == HAP_LABLE)
+				{/****************************may have hap bugs********************************/
 					av[i].del = 1;
 					asg_arc_del(g, av[i].v^1, v^1, 1);
 				}
@@ -313,6 +321,36 @@ typedef struct {
 	uint32_t pre_n_seq, seqID; 
 } C_graph;
 
+typedef struct {
+	kvec_t(uint32_t) a;
+	uint32_t i;
+} kvec_t_u32_warp;
+
+typedef struct {
+	kvec_t(uint64_t) a;
+	uint64_t i;
+} kvec_t_u64_warp;
+
+
+typedef struct {
+	uint32_t q_pos;
+	uint32_t t_pos;
+	uint32_t t_id;
+	uint32_t is_color;
+} Hap_Align;
+
+typedef struct {
+	kvec_t(Hap_Align) x;
+	uint64_t i;
+} Hap_Align_warp;
+
+typedef struct {
+	buf_t* b_0;
+    uint32_t untigI;
+    uint32_t readI;
+    uint32_t offset;
+} rIdContig;
+
 // count the number of outgoing arcs, including reduced arcs
 static inline int count_out_with_del(const asg_t *g, uint32_t v)
 {
@@ -350,4 +388,38 @@ void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long lo
 void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, 
 uint64_t* source_index, long long listLen);
 void print_revise_edges(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen);
+
+#define EvaluateLen(U, id) ((U).a[(id)].start)
+#define IsMerge(U, id) ((U).a[(id)].end)
+#define kv_reuse(v, rn, rm, r) ((v).n = (rn), (v).m = (rm), (v).a = (r))
+#define long_tip(U, id, threshold) ((EvaluateLen((U), (id))>=(threshold))&&(!((U).a[(id)].circ)))
+///there are threee cases: 
+///1. if this untig is too long (>maxShortUntig), it must be not short untig/must be a long untig
+///2. if this untig is long (>minLongUntig && EvaluateLen(ug->u, av[i].v>>1) > (EvaluateLen(ug->u, v>>1)*l_untig_rate)), it might be a long tip
+#define check_long_tip(U, id, minLongUntig, maxShortUntig, ShortUntigRate, mainLen) \
+                        ((!((U).a[(id)].circ)) \
+                        && \
+                        ((EvaluateLen((U), (id)) > (maxShortUntig))\
+                        ||\
+                        ((long_tip((U), (id), (minLongUntig)))\
+                        &&\
+                        (EvaluateLen((U), (id)) > (ShortUntigRate)*(mainLen)))))
+#define Get_vis(visit, v, d) (((visit)[(v)>>1])&(((((v)<<(d))&1)+1)))
+#define Set_vis(visit, v, d) (((visit)[(v)>>1])|=(((((v)<<(d))&1)+1)))
+
+
+
+
+
+typedef struct {
+	uint64_t len;
+	uint32_t* index;
+} R_to_U;
+
+void init_R_to_U(R_to_U* x, uint64_t len);
+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);
+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);
+
 #endif

README.md

@@ -43,9 +43,9 @@ assembly in a few hours. Hifiasm has been tested on the following datasets:
 
 |<sub>Dataset<sub>|<sub>GSize<sub>|<sub>Cov<sub>|<sub>Asm options<sub>|<sub>CPU time<sub>|<sub>Wall time<sub>|<sub>RAM<sub>|<sub>[unitig][unitig]/[contig][unitig] N50<sup>[1]</sup><sub>|
 |:---------------|-----:|-----:|:---------------------|-------:|--------:|----:|----------------:|
-|<sub>[Human NA12878]<sub>|<sub>3Gb<sub>|<sub>x28<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>200h<sub>|    <sub>5h32m<sub>|<sub>114G<sub>|<sub>93.5Kb/21.5Mb<sub>|
+|<sub>[Human NA12878]<sub>|<sub>3Gb<sub>|<sub>x28<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>200h<sub>|    <sub>5h32m<sub>|<sub>114G<sub>|<sub>93.5Kb/28.2Mb<sub>|
-|<sub>[Human HG002]<sub>|<sub>3Gb<sub>|<sub>x43<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>405h10m<sub>|<sub>12h7m<sub>|<sub>146G<sub>|<sub>320kb/31.9Mb<sub>|
+|<sub>[Human HG002]<sub>|<sub>3Gb<sub>|<sub>x43<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>405h10m<sub>|<sub>12h7m<sub>|<sub>146G<sub>|<sub>320kb/35.3Mb<sub>|
-|<sub>[Human CHM13]<sub>|<sub>3Gb<sub>|<sub>x27<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>157h28m<sub>|<sub>5h10m<sub>|<sub>85.8G<sub>|<sub>NA<sup>[2]</sup>/39.8Mb<sub>|
+|<sub>[Human CHM13]<sub>|<sub>3Gb<sub>|<sub>x27<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>157h28m<sub>|<sub>5h10m<sub>|<sub>85.8G<sub>|<sub>NA<sup>[2]</sup>/41.4Mb<sub>|
 |<sub>[Butterfly]<sub>|<sub>358Mb<sub>|<sub>x35<sub>|<sub>-k 40 -t 42 -r 2 -z 20<sub>|<sub>17h6m<sub>|<sub>36m<sub>|<sub>16G<sub>|<sub>7.5Mb/NA<sup>[3]</sup><sub>|
 
 <sub>[1] unitig N50 is the N50 of assembly graph with haplotype information (i.e., bubbles), while the contig N50 is the N50 of haplotype collapsed assembly (i.e., without bubbles).