xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-10-07 02:08:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into dev-lh3

Browse Source
This commit is contained in:
Heng Li
2021-03-18 16:11:43 -04:00
parent 8e98fb2276 e6e6dbf7b3
commit 863b20773a
8 changed files with 5650 additions and 2397 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+24 -13
View File
@@ -23,7 +23,7 @@ static ko_longopt_t long_options[] = {
{ "ex-iter", ko_required_argument, 308 },
{ "purge-cov", ko_required_argument, 309 },
{ "pri-range", ko_required_argument, 310 },
{ "high-het", ko_no_argument, 311 },
///{ "high-het", ko_no_argument, 311 },
{ "lowQ", ko_required_argument, 312 },
{ "min-hist-cnt", ko_required_argument, 313 },
{ "h1", ko_required_argument, 314 },
@@ -89,14 +89,14 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " -4 FILE list of hap2/maternal read names []\n");
fprintf(stderr, " Purge-dups:\n");
fprintf(stderr, " -l INT purge level. 0: no purging; 1: light; 2: aggressive [0 for trio; 2 for unzip]\n");
fprintf(stderr, " -s FLOAT similarity threshold for duplicate haplotigs [%g]\n",
asm_opt->purge_simi_rate);
fprintf(stderr, " -l INT purge level. 0: no purging; 1: light; 2/3: aggressive [0 for trio; 2 for unzip]\n");
fprintf(stderr, " -s FLOAT similarity threshold for duplicate haplotigs [%g for -l1/-l2, %g for -l3]\n",
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, " coverage upper bound of Purge-dups [auto]\n");
fprintf(stderr, " --high-het enable this mode for high heterozygosity sample [experimental, not stable]\n");
///fprintf(stderr, " --high-het enable this mode for high heterozygosity sample [experimental, not stable]\n");
fprintf(stderr, " Hi-C-partition [experimental, not stable]:\n");
fprintf(stderr, " --h1 FILEs file names of Hi-C R1 [r1_1.fq,r1_2.fq,...]\n");
@@ -151,10 +151,11 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->mid_cnt = 5;
asm_opt->purge_level_primary = 2;
asm_opt->purge_level_trio = 0;
asm_opt->purge_simi_rate = 0.75;
asm_opt->purge_simi_rate_hic = 0.85;
asm_opt->purge_simi_rate_l2 = 0.75;
asm_opt->purge_simi_rate_l3 = 0.55;
///asm_opt->purge_simi_rate_hic = 0.85;
asm_opt->purge_overlap_len = 1;
asm_opt->purge_overlap_len_hic = 50;
///asm_opt->purge_overlap_len_hic = 50;
asm_opt->recover_atg_cov_min = -1024;
asm_opt->recover_atg_cov_max = INT_MAX;
asm_opt->hom_global_coverage = -1;
@@ -165,6 +166,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->b_low_cov = 0;
asm_opt->b_high_cov = -1;
asm_opt->m_rate = 0.75;
asm_opt->hap_occ = 1;
}
void destory_enzyme(enzyme* f)
@@ -361,9 +363,9 @@ int check_option(hifiasm_opt_t* asm_opt)
return 0;
}
if(asm_opt->purge_level_primary < 0 || asm_opt->purge_level_primary > 2)
if(asm_opt->purge_level_primary < 0 || asm_opt->purge_level_primary > 3)
{
fprintf(stderr, "[ERROR] the level of purge-dup should be [0, 2] (-l)\n");
fprintf(stderr, "[ERROR] the level of purge-dup should be [0, 3] (-l)\n");
return 0;
}
@@ -454,6 +456,12 @@ int check_option(hifiasm_opt_t* asm_opt)
return 0;
}
if(asm_opt->purge_simi_thres < 0)
{
fprintf(stderr, "[ERROR] [-s] must >= 0\n");
return 0;
}
// fprintf(stderr, "input file num: %d\n", asm_opt->num_reads);
// fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
// fprintf(stderr, "number of threads: %d\n", asm_opt->thread_num);
@@ -471,7 +479,7 @@ int check_option(hifiasm_opt_t* asm_opt)
// fprintf(stderr, "mid_cnt: %d\n", asm_opt->mid_cnt);
// fprintf(stderr, "purge_level_primary: %d\n", asm_opt->purge_level_primary);
// fprintf(stderr, "purge_level_trio: %d\n", asm_opt->purge_level_trio);
// fprintf(stderr, "purge_simi_rate: %f\n", asm_opt->purge_simi_rate);
// fprintf(stderr, "purge_simi_thres: %f\n", asm_opt->purge_simi_thres);
// fprintf(stderr, "purge_overlap_len: %d\n", asm_opt->purge_overlap_len);
return 1;
@@ -621,7 +629,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
asm_opt->recover_atg_cov_min = asm_opt->recover_atg_cov_max = -1;
}
}
else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET;
///else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET;
else if (c == 312) asm_opt->bed_inconsist_rate = atoi(opt.arg);
else if (c == 313) asm_opt->min_hist_kmer_cnt = atoi(opt.arg);
else if (c == 314) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[0]), 0);
@@ -634,7 +642,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
{ ///0: disable purge_dup; 1: purge containment; 2: purge overlap
asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);
}
else if (c == 's') asm_opt->purge_simi_rate = atof(opt.arg);
else if (c == 's') asm_opt->purge_simi_rate_l2 = asm_opt->purge_simi_rate_l3 = atof(opt.arg);
else if (c == 'O') asm_opt->purge_overlap_len = atoll(opt.arg);
else if (c == ':')
{
@@ -648,6 +656,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
}
}
if(asm_opt->purge_level_primary > 2) asm_opt->purge_simi_thres = asm_opt->purge_simi_rate_l3;
else asm_opt->purge_simi_thres = asm_opt->purge_simi_rate_l2;
if (argc == opt.ind)
{
CommandLines.h
+7 -4
View File
@@ -3,7 +3,7 @@
#include <pthread.h>
#define HA_VERSION "0.14-r312"
#define HA_VERSION "0.14.1-r314"
#define VERBOSE 0
@@ -72,7 +72,7 @@ typedef struct {
int purge_level_primary;
int purge_level_trio;
int purge_overlap_len;
int purge_overlap_len_hic;
///int purge_overlap_len_hic;
int recover_atg_cov_min;
int recover_atg_cov_max;
int hom_global_coverage;
@@ -82,8 +82,10 @@ typedef struct {
float max_hang_rate;
float min_drop_rate;
float max_drop_rate;
float purge_simi_rate;
float purge_simi_rate_hic;
float purge_simi_rate_l2;
float purge_simi_rate_l3;
float purge_simi_thres;
///float purge_simi_rate_hic;
long long small_pop_bubble_size;
long long large_pop_bubble_size;
@@ -92,6 +94,7 @@ typedef struct {
long long num_recorrected_bases;
long long mem_buf;
long long coverage;
int hap_occ;
} hifiasm_opt_t;
Overlaps.cpp
+1922 -1037
View File
File diff suppressed because it is too large Load Diff
Overlaps.h
+84 -73
View File
@@ -478,7 +478,6 @@ void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8
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);
void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, kvec_asg_arc_t_warp* edge, int require_equal_nv, int test_tangle);
int asg_pop_bubble_primary(asg_t *g, int max_dist);
long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *g, long long circleLen, int is_drop);
typedef struct {
@@ -493,9 +492,21 @@ typedef struct {
uint32_t new_edges_i;
} Edge_iter;
typedef struct {
asg_arc_t x;
uint64_t Off;
uint64_t weight;
}asg_arc_t_offset;
typedef struct {
kvec_t(asg_arc_t_offset) a;
uint64_t i;
}kvec_asg_arc_t_offset;
void init_Edge_iter(asg_t* g, uint32_t v, asg_arc_t* new_edges, uint32_t new_edges_n, Edge_iter* x);
int get_arc_t(Edge_iter* x, asg_arc_t* get);
int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag);
inline int get_real_length(asg_t *g, uint32_t v, uint32_t* v_s)
@@ -540,55 +551,6 @@ inline uint32_t check_tip(asg_t *sg, uint32_t begNode, uint32_t* endNode, buf_t*
}
}
inline uint32_t get_unitig_back(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode,
long long* nodeLen, long long* baseLen, buf_t* b)
{
ma_utg_v* u = NULL;
uint32_t v = begNode, w, k;
uint32_t kv;
(*nodeLen) = (*baseLen) = 0;
(*endNode) = (uint32_t)-1;
if(ug!=NULL) u = &(ug->u);
while (1)
{
kv = get_real_length(sg, v, NULL);
(*endNode) = v;
if(u == NULL)
{
(*nodeLen)++;
}
else
{
(*nodeLen) += EvaluateLen((*u), v>>1);
}
if(b) kv_push(uint32_t, b->b, v);
///means reach the end of a unitig
if(kv!=1) (*baseLen) += sg->seq[v>>1].len;
if(kv==0) return END_TIPS;
if(kv>1) return MUL_OUTPUT;
///kv must be 1 here
kv = get_real_length(sg, v, &w);
///means reach the end of a unitig
if(get_real_length(sg, w^1, NULL)!=1)
{
(*baseLen) += sg->seq[v>>1].len;
return MUL_INPUT;
}
for (k = 0; k < asg_arc_n(sg, v); k++)
{
if(asg_arc_a(sg, v)[k].del) continue;
///here is just one undeleted edge
(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
break;
}
v = w;
if(v == begNode) return LOOP;
}
}
inline uint32_t get_unitig(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode,
long long* nodeLen, long long* baseLen, long long* max_stop_nodeLen, long long* max_stop_baseLen,
uint32_t stops_threshold, buf_t* b)
@@ -1019,31 +981,49 @@ typedef struct {
uint32_t total;
} Trio_counter;
typedef struct {
uint32_t p; // the optimal parent vertex
uint32_t d; // the shortest distance from the initial vertex
uint32_t r:31, s:1; // r: the number of remaining incoming arc; s: state
} binfo_s_t;
typedef struct {
///all information for each node
binfo_s_t *a;
kvec_t(uint32_t) S; // set of vertices without parents, nodes with all incoming edges visited
kvec_t(uint32_t) b; // visited vertices
kvec_t(uint32_t) e; // visited edges/arcs
} buf_s_t;
typedef struct{
buf_s_t *b;
uint32_t n_thres, n_reads;
asg_t *g;
uint32_t check_cross;
uint64_t bub_dist;
} bub_label_t;
void resolve_tangles(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, uint32_t trio_flag,
float drop_ratio);
void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
void rescue_contained_reads_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t chainLenThres, uint32_t is_bubble_check,
uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes);
R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t chainLenThres, uint32_t is_bubble_check,
uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes, bub_label_t* b_mask_t);
void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
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);
R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t is_bubble_check, uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t);
void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate);
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);
void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges);
ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges, bub_label_t* b_mask_t);
void get_unitig_trio_flag(ma_utg_t* nsu, uint32_t flag, uint32_t* require, uint32_t* non_require, uint32_t* ambigious);
void rescue_missing_overlaps_backward(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t backward_steps,
uint32_t is_bubble_check, uint32_t is_primary_check);
R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t backward_steps, uint32_t is_bubble_check, uint32_t is_primary_check, bub_label_t* b_mask_t);
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, uint64_t* path_base_len, uint64_t* path_nodes);
int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio);
void asg_bub_backtrack_primary(asg_t *g, uint32_t v0, buf_t *b);
typedef struct{
double weight;
@@ -1072,10 +1052,38 @@ typedef struct{
kvec_t(hc_linkeage) a;
kvec_t(uint64_t) enzymes;
kvec_t(bed_in) bed;
uint32_t* u_idx;
uint64_t r_num;
///uint32_t* u_idx;
///uint64_t r_num;
} hc_links;
typedef struct{
kvec_t(uint32_t) uIDs;
kvec_t(uint32_t) iDXs;
kvec_t(uint32_t) rescue_hom;
uint32_t* u_idx;
uint32_t r_num;
uint32_t chain_num;
uint32_t l0_chain, l1_chain;
}trans_chain;
typedef struct {
uint32_t n;
uint32_t* cov;
uint64_t* pos_idx;
ma_hit_t_alloc* reverse_sources;
ma_sub_t *coverage_cut;
R_to_U* ruIndex;
asg_t *read_g;
int max_hang;
int min_ovlp;
kvec_asg_arc_t_offset u_buffer;
kvec_t_i32_warp tailIndex;
kvec_t_i32_warp prevIndex;
///hc_links* link;
trans_chain* t_ch;
}hap_cov_t;
typedef struct{
///kvec_t(hc_edge) a;
size_t n, m;
@@ -1084,22 +1092,22 @@ typedef struct{
void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num);
void destory_hc_links(hc_links* link);
void clean_primary_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold,
R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen,
uint32_t miniBiGraph, float chimeric_rate, int is_final_clean, int just_bubble_pop,
float drop_ratio, hc_links* link);
uint64_t get_bub_pop_max_dist(asg_t *g, buf_t *b);
uint64_t get_bub_pop_max_dist_advance(asg_t *g, buf_t *b);
int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov);
uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, uint64_t max_dist, buf_t *b, uint32_t positive_flag,
uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes, hap_cov_t *cov);
void adjust_utg_by_primary(ma_ug_t **ug, asg_t* read_g, float drop_rate,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut,
long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold,
long long tipsLen, float tip_drop_ratio, long long stops_threshold,
R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link);
void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg);
kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link, bub_label_t* b_mask_t);
ma_ug_t* copy_untig_graph(ma_ug_t *src);
ma_ug_t* output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name,
uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist,
uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex,
float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench);
float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench, bub_label_t* b_mask_t);
asg_t* copy_read_graph(asg_t *src);
ma_ug_t *ma_ug_gen(asg_t *g);
void ma_ug_destroy(ma_ug_t *ug);
@@ -1112,6 +1120,9 @@ inline int inter_interval(int a_s, int a_e, int b_s, int b_e, int* i_s, int* i_e
return 1;
}
#define JUNK_COV 5
#define DISCARD_RATE 0.8
Purge_Dups.cpp
+1391 -1049
View File
File diff suppressed because it is too large Load Diff
Purge_Dups.h
+8 -2
View File
@@ -11,14 +11,20 @@
#define HET_PEAK_RATE (HOM_PEAK_RATE*2)
#define ALTER_COV_THRES 0.9
#define REAL_ALTER_THRES 0.1
#define CHAIN_FILTER_RATE 0.7
void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density,
uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio,
uint32_t just_contain, uint32_t just_coverage, hc_links* link);
uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain,
uint32_t just_coverage, hap_cov_t *cov);
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 get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen);
void enable_debug_mode(uint32_t mode);
hap_cov_t* init_hap_cov_t(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* sources, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources,
ma_sub_t *coverage_cut, int max_hang, int min_ovlp, uint32_t is_collect_trans);
void destory_hap_cov_t(hap_cov_t **x);
void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd);
#endif
hic.cpp
+2210 -212
View File
File diff suppressed because it is too large Load Diff
hifiasm.1
+4 -7
View File
@@ -286,12 +286,14 @@ times in the other sample.
.TP 10
.BI -l \ INT
Level of purge-dup. 0 to disable purge-dup, 1 to only purge contained haplotigs,
2 to purge all types of haplotigs. In default, [2] for non-trio assembly, [0] for trio assembly.
2 to purge all types of haplotigs, 3 to purge all types of haplotigs in most aggressive way
for high heterozygosity sample.
In default, [2] for non-trio assembly, [0] for trio assembly.
For trio assembly, only level 0 and level 1 are allowed.
.TP
.BI -s \ FLOAT
Similarity threshold for duplicate haplotigs that should be purged [0.75].
Similarity threshold for duplicate haplotigs that should be purged [0.75 for -l1/-l2, 0.55 for -l3].
.TP
.BI -O \ FLOAT
@@ -302,11 +304,6 @@ 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