Merge branch 'master' of https://github.com/chhylp123/Long_read_assembly

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;
@@ -46,6 +46,7 @@ void Print_H(hifiasm_opt_t* asm_opt)
 
 void init_opt(hifiasm_opt_t* asm_opt)
 {
+    asm_opt->coverage = -1;
     asm_opt->num_reads = 0;
     asm_opt->read_file_names = NULL;
     asm_opt->output_file_name = (char*)(DEFAULT_OUTPUT);

CommandLines.h

@@ -36,6 +36,7 @@ typedef struct {
     long long num_bases;
     long long num_corrected_bases;
     long long num_recorrected_bases;
+    long long coverage;
 } hifiasm_opt_t;
 
 extern 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;
@@ -102,6 +108,8 @@ typedef struct {
 	uint32_t m_arc, n_arc:31, is_srt:1;
 	asg_arc_t *arc;
 	uint32_t m_seq, n_seq:31, is_symm:1;
+	uint32_t r_seq;
+
 	asg_seq_t *seq;
 	uint64_t *idx;
 
@@ -239,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) 
@@ -252,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);
 				}
@@ -304,6 +314,43 @@ typedef struct {
 	kvec_t(uint32_t) e; // visited edges/arcs
 } buf_t;
 
+
+typedef struct {
+	kvec_t(uint64_t) Nodes; 
+	kvec_t(uint64_t) Edges; 
+	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)
 {
@@ -341,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

@@ -10,7 +10,7 @@ cd hifiasm && make
 
 ## Introduction
 
-Hifiasm is a fast haplotype-reserved de novo assembler for PacBio
+Hifiasm is a fast haplotype-resolved de novo assembler for PacBio
 Hifi reads. Unlike most existing assemblers, hifiasm starts from uncollapsed
 genome. Thus, it is able to keep the haplotype information as much as possible.
 The input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and its
@@ -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/18.6Mb<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/29.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 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/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>/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).