seperate purgeing and hom peak

CommandLines.cpp

@@ -21,7 +21,7 @@ static ko_longopt_t long_options[] = {
 	{ "max-od-final",  ko_no_argument, 306 },
 	{ "ex-list",       ko_required_argument, 307 },
 	{ "ex-iter",       ko_required_argument, 308 },
-    { "purge-cov",     ko_required_argument, 309 },
+    { "hom-cov",     ko_required_argument, 309 },
     { "pri-range",     ko_required_argument, 310 },
     { "lowQ",          ko_required_argument, 312 },
 	{ "min-hist-cnt",  ko_required_argument, 313 },
@@ -39,6 +39,7 @@ static ko_longopt_t long_options[] = {
     { "n-hap",      ko_required_argument, 325 },
     { "n-weight",      ko_required_argument, 326 },
     { "l-msjoin",      ko_required_argument, 327 },
+    { "purge-max",     ko_required_argument, 328 },
 	{ 0, 0, 0 }
 };
 
@@ -55,9 +56,8 @@ void Print_H(hifiasm_opt_t* asm_opt)
     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, "    -h           show help information\n");
     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);
@@ -68,12 +68,16 @@ void Print_H(hifiasm_opt_t* asm_opt)
     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, "    -z INT       length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
     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, "    -i           ignore saved read correction and overlaps\n");
     fprintf(stderr, "    -u           disable post join contigs step which may improve N50\n");
+    fprintf(stderr, "    --hom-cov    INT\n");
+    fprintf(stderr, "                 homozygous read coverage [auto]\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");
@@ -105,7 +109,7 @@ void Print_H(hifiasm_opt_t* asm_opt)
                                       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, "    --purge-max  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);
@@ -183,6 +187,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->recover_atg_cov_max = INT_MAX;
     asm_opt->hom_global_coverage = -1;
     asm_opt->hom_global_coverage_set = 0;
+    asm_opt->pur_global_coverage = -1;
     asm_opt->bed_inconsist_rate = 70;
     asm_opt->hic_inconsist_rate = 30;
     ///asm_opt->bub_mer_length = 3;
@@ -410,7 +415,13 @@ int check_option(hifiasm_opt_t* asm_opt)
 
     if(asm_opt->hom_global_coverage < 0 && asm_opt->hom_global_coverage != -1)
     {
-        fprintf(stderr, "[ERROR] purge duplication coverage threshold should be >= 0 (--purge-cov)\n");
+        fprintf(stderr, "[ERROR] homozygous read coverage should be >= 0 (--hom-cov)\n");
+        return 0;
+    }
+
+    if(asm_opt->pur_global_coverage < 0 && asm_opt->pur_global_coverage != -1)
+    {
+        fprintf(stderr, "[ERROR] purge duplication coverage threshold should be >= 0 (--purge-max)\n");
         return 0;
     }
 
@@ -663,6 +674,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 325) asm_opt->polyploidy = atoi(opt.arg);
         else if (c == 326) asm_opt->n_weight = atoi(opt.arg);
         else if (c == 327) asm_opt->misjoin_len = atol(opt.arg);
+        else if (c == 328) asm_opt->pur_global_coverage = 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);

CommandLines.h

@@ -4,7 +4,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.15.4-r348"
+#define HA_VERSION "0.15.4-r349"
 
 #define VERBOSE 0
 
@@ -77,6 +77,7 @@ typedef struct {
     int recover_atg_cov_max;
     int hom_global_coverage;
     int hom_global_coverage_set;
+    int pur_global_coverage;
     int bed_inconsist_rate;
     int hic_inconsist_rate;
 

Overlaps.cpp

@@ -16995,7 +16995,8 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
 
     if(asm_opt.recover_atg_cov_min == -1024)
     {
-        asm_opt.recover_atg_cov_max = asm_opt.hom_global_coverage/HOM_PEAK_RATE;
+        asm_opt.recover_atg_cov_max = (asm_opt.hom_global_coverage_set?
+            (asm_opt.hom_global_coverage):(((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE)));         
         asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.85;
         asm_opt.recover_atg_cov_max = INT32_MAX;
     }
@@ -24806,7 +24807,8 @@ uint32_t collect_p_trans, uint32_t collect_p_trans_f)
 
     if(asm_opt.recover_atg_cov_min == -1024)
     {
-        asm_opt.recover_atg_cov_max = asm_opt.hom_global_coverage/HOM_PEAK_RATE;
+        asm_opt.recover_atg_cov_max = (asm_opt.hom_global_coverage_set?
+            (asm_opt.hom_global_coverage):(((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE)));
         asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.85;
         asm_opt.recover_atg_cov_max = INT32_MAX;
     }

Purge_Dups.cpp

@@ -5267,9 +5267,14 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle
     hap_alignment_struct_pip hap_buf;
     long long k_mer_only, coverage_only;
 
-    if(asm_opt.hom_global_coverage != -1)
+    if(asm_opt.pur_global_coverage != -1)
     {
-        hap_buf.cov_threshold = asm_opt.hom_global_coverage;
+        hap_buf.cov_threshold = asm_opt.pur_global_coverage;
+    }
+    else if(asm_opt.hom_global_coverage != -1)
+    {
+        hap_buf.cov_threshold = (asm_opt.hom_global_coverage_set?
+                (((double)asm_opt.hom_global_coverage)*((double)HOM_PEAK_RATE)):(asm_opt.hom_global_coverage));
     }
     else
     {

hic.cpp

@@ -43,6 +43,7 @@ KRADIX_SORT_INIT(u_trans_m, u_trans_t, u_trans_m_key, 8)
 KRADIX_SORT_INIT(u_trans_occ, u_trans_t, u_trans_occ_key, member_size(u_trans_t, occ))
 
 #define is_hom_hit(a) ((a).id == (uint64_t)-1)
+#define HC_PT_MA 65
 
 typedef struct {
     kv_gg_status sg;
@@ -555,7 +556,9 @@ int write_hc_pt_index(ha_ug_index* idx, char* file_name)
         free(gfa_name);
         return 0;
     }
+    uint64_t i = HC_PT_MA;
 
+    fwrite(&i, sizeof(i), 1, fp);
     fwrite(&idx->uID_bits, sizeof(idx->uID_bits), 1, fp);
     fwrite(&idx->uID_mode, sizeof(idx->uID_mode), 1, fp);
     fwrite(&idx->pos_bits, sizeof(idx->pos_bits), 1, fp);
@@ -566,7 +569,6 @@ int write_hc_pt_index(ha_ug_index* idx, char* file_name)
     fwrite(&idx->tot, sizeof(idx->tot), 1, fp);
     fwrite(&idx->tot_pos, sizeof(idx->tot_pos), 1, fp);
 
-    uint64_t i = 0;
     for (i = 0; i < idx->tot; i++)
     {
         fwrite(&idx->idx_buf[i].n, sizeof(idx->idx_buf[i].n), 1, fp);
@@ -596,7 +598,19 @@ int load_hc_pt_index(ha_ug_index** r_idx, ma_ug_t *ug, char* file_name)
         return 0;
     }
     ha_ug_index* idx = NULL; CALLOC(idx, 1);
+    uint64_t i;
 
+    flag += fread(&i, sizeof(i), 1, fp);
+    if(i != HC_PT_MA)
+    {
+        free(gfa_name);
+        destory_hc_pt_index(idx);
+        free(idx);
+        (*r_idx) = NULL;
+        fclose(fp);
+        fprintf(stderr, "[M::%s::] ==> Renew Hi-C index\n", __func__);
+        return 0;
+    }
     flag += fread(&idx->uID_bits, sizeof(idx->uID_bits), 1, fp);
     flag += fread(&idx->uID_mode, sizeof(idx->uID_mode), 1, fp);
     flag += fread(&idx->pos_bits, sizeof(idx->pos_bits), 1, fp);
@@ -607,7 +621,7 @@ int load_hc_pt_index(ha_ug_index** r_idx, ma_ug_t *ug, char* file_name)
     flag += fread(&idx->tot, sizeof(idx->tot), 1, fp);
     flag += fread(&idx->tot_pos, sizeof(idx->tot_pos), 1, fp);
     MALLOC(idx->idx_buf, idx->tot);
-    uint64_t i = 0;
+
     for (i = 0; i < idx->tot; i++)
     {
         flag += fread(&idx->idx_buf[i].n, sizeof(idx->idx_buf[i].n), 1, fp);
@@ -877,6 +891,14 @@ ha_ug_index* build_unitig_index(ma_ug_t *ug, int k, uint64_t up_occ, uint64_t lo
 
     fprintf(stderr, "[M::%s::%.3f] ==> HiC index has been built\n", __func__, yak_realtime()-index_time);
     
+    uint64_t i;
+    for (i = 0; i < idx->tot; i++)
+    {
+        kv_destroy(pl.cnt[i].a);
+        kv_destroy(pl.buf[i].a);
+    }
+    free(pl.cnt); free(pl.buf);
+
     return idx;
 }
 
@@ -4043,6 +4065,7 @@ void update_containment_distance(asg_t *sg, kv_u_trans_t *ta, hc_links* link)
     }
 
     kv_destroy(buf.a);
+    free(uc_idx);
 }
 
 
@@ -4877,7 +4900,12 @@ int load_hc_hits(kvec_pe_hit* hits, ma_ug_t* ug, const char *fn)
 
     FILE* fp = NULL; 
     fp = fopen(buf, "r"); 
-    if(!fp) return 0;
+    if(!fp) 
+    {
+        free(buf);
+        return 0;
+    }
+    
 
     kv_init(hits->a);
     flag += fread(&hits->a.n, sizeof(hits->a.n), 1, fp);
@@ -9511,6 +9539,15 @@ H_partition* hap, int8_t *s, mc_gg_status *sa, trans_idx* dis)
     trans_p_t* p = NULL;
     dis->n = 0;
     uint64_t bin_size = MIN(2250, buf.n>>8), m;
+    if(bin_size == 0)
+    {
+        for (i = 8; i > 0; i--)
+        {
+            bin_size = buf.n>>i;
+            if(bin_size > 0) break;
+        }
+        if(bin_size == 0) bin_size = buf.n;
+    }
     i = 0;
     while (i < buf.n)
     {
@@ -16448,7 +16485,7 @@ void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt,
     else hic_short_align_poy(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt);
     
     
-    destory_hc_pt_index(ug_index);
+    destory_hc_pt_index(ug_index);free(ug_index);
 }
 
 

khashl.h

@@ -136,25 +136,33 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 #define __KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \
 	SCOPE khint_t prefix##_save(HType *h, FILE* fp) { \
 		if (!h) return 0; \
+		uint8_t ff; \
 		khint_t n_buckets = (h->keys? 1U<<h->bits : 0U); \
 		fwrite(&n_buckets, sizeof(n_buckets), 1, fp); \
 		fwrite(&h->bits, sizeof(h->bits), 1, fp); \
 		fwrite(&h->count, sizeof(h->count), 1, fp); \
-		fwrite(h->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
+		ff = h->used? 1:0; fwrite(&ff, sizeof(ff), 1, fp); \
-		fwrite(h->keys, sizeof(khkey_t), n_buckets, fp); \
+		if(ff) fwrite(h->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
+		ff = h->keys? 1:0; fwrite(&ff, sizeof(ff), 1, fp); \
+		if(ff) fwrite(h->keys, sizeof(khkey_t), n_buckets, fp); \
 		return 1; \
 	} \
 	SCOPE khint_t prefix##_load(HType **h, FILE* fp) { \
 		(*h) = prefix##_init(); \
 		khint_t n_buckets; \
 		uint64_t flag = 0;\
+		uint8_t ff; \
 		flag += fread(&n_buckets, sizeof(n_buckets), 1, fp); \
 		flag += fread(&(*h)->bits, sizeof((*h)->bits), 1, fp); \
 		flag += fread(&(*h)->count, sizeof((*h)->count), 1, fp); \
+		flag += fread(&ff, sizeof(ff), 1, fp); \
+		if(ff) {\
 		(*h)->used = (khint32_t*)kmalloc(__kh_fsize(n_buckets) * sizeof(khint32_t)); \
+		flag += fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); }\
+		flag += fread(&ff, sizeof(ff), 1, fp); \
+		if(ff) {\
 		(*h)->keys = (khkey_t*)kmalloc(n_buckets * sizeof(khkey_t)); \
-		flag += fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
+		flag += fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); }\
-		flag += fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \
 		return 1; \
 	} \