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

for v0.14

Browse Source
This commit is contained in:
chhylp123
2021-02-09 13:00:23 -05:00
parent fda13cb0e4
commit 8467d58ef1
7 changed files with 2928 additions and 925 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+95 -45
View File
@@ -29,6 +29,7 @@ static ko_longopt_t long_options[] = {
{ "h1", ko_required_argument, 314 },
{ "h2", ko_required_argument, 315 },
{ "enzyme", ko_required_argument, 316 },
{ "b-cov", ko_required_argument, 317 },
{ 0, 0, 0 }
};
@@ -66,6 +67,9 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " -u disable post join contigs step which may improve N50\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");
fprintf(stderr, " break contigs at breakpoints with coverage drop at <INT-fold coverage [%d]\n", asm_opt->break_cov);
// fprintf(stderr, " --pri-range INT1[,INT2]\n");
// fprintf(stderr, " keep contigs with coverage in this range in p_ctg.gfa; -1 to disable [auto,inf]\n");
@@ -87,6 +91,9 @@ void Print_H(hifiasm_opt_t* asm_opt)
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, " Hi-C-partition [experimental, not stable]:\n");
fprintf(stderr, " --h1 FILEs file names of Hi-C R1 [r1_1.fq,r1_2.fq,...]\n");
fprintf(stderr, " --h2 FILEs file names of Hi-C R2 [r2_1.fq,r2_2.fq,...]\n");
fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
fprintf(stderr, "See `man ./hifiasm.1' for detailed description of these command-line options.\n");
@@ -101,9 +108,9 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->read_file_names = NULL;
asm_opt->output_file_name = (char*)(DEFAULT_OUTPUT);
asm_opt->required_read_name = NULL;
asm_opt->hic_files[0] = NULL;
asm_opt->hic_files[1] = NULL;
asm_opt->hic_enzymes = NULL;
asm_opt->hic_reads[0] = NULL;
asm_opt->hic_reads[1] = NULL;
asm_opt->thread_num = 1;
asm_opt->k_mer_length = 51;
asm_opt->hic_mer_length = 31;
@@ -146,22 +153,30 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->hic_inconsist_rate = 30;
///asm_opt->bub_mer_length = 3;
asm_opt->bub_mer_length = 1000000;
asm_opt->break_cov = 0;
}
void destory_enzyme(enzyme* f)
{
int i;
if(f != NULL)
{
for (i = 0; i < f->n; i++)
{
free(f->a[i]);
}
free(f->a);
free(f->l);
free(f);
}
}
void destory_opt(hifiasm_opt_t* asm_opt)
{
if(asm_opt->read_file_names != NULL) free(asm_opt->read_file_names);
if(asm_opt->hic_enzymes != NULL)
{
int i;
for (i = 0; i < asm_opt->hic_enzymes->n; i++)
{
free(asm_opt->hic_enzymes->a[i]);
}
free(asm_opt->hic_enzymes->a);
free(asm_opt->hic_enzymes->l);
free(asm_opt->hic_enzymes);
}
if(asm_opt->hic_enzymes != NULL) destory_enzyme(asm_opt->hic_enzymes);
if(asm_opt->hic_reads[0] != NULL) destory_enzyme(asm_opt->hic_reads[0]);
if(asm_opt->hic_reads[1] != NULL) destory_enzyme(asm_opt->hic_reads[1]);
}
void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round)
@@ -201,6 +216,16 @@ static int check_file(char* name, const char* opt)
return 1;
}
static int check_hic_reads(enzyme* f, const char* opt)
{
int i;
for (i = 0; i < f->n; i++)
{
if(check_file(f->a[i], opt) == 0) return 0;
}
return 1;
}
int check_option(hifiasm_opt_t* asm_opt)
{
if(asm_opt->read_file_names == NULL || asm_opt->num_reads == 0)
@@ -355,25 +380,45 @@ int check_option(hifiasm_opt_t* asm_opt)
if(asm_opt->fn_bin_list[0] != NULL && check_file(asm_opt->fn_bin_list[0], "LIST1") == 0) return 0;
if(asm_opt->fn_bin_list[1] != NULL && check_file(asm_opt->fn_bin_list[1], "LIST2") == 0) return 0;
if(asm_opt->required_read_name != NULL && check_file(asm_opt->required_read_name, "b") == 0) return 0;
if(asm_opt->hic_files[0] != NULL && check_file(asm_opt->hic_files[0], "HIC1") == 0) return 0;
if(asm_opt->hic_files[1] != NULL && check_file(asm_opt->hic_files[1], "HIC2") == 0) return 0;
if(asm_opt->hic_files[0] != NULL && asm_opt->hic_files[1] == NULL)
if(asm_opt->hic_reads[0] != NULL && check_hic_reads(asm_opt->hic_reads[0], "HIC1") == 0) return 0;
if(asm_opt->hic_reads[1] != NULL && check_hic_reads(asm_opt->hic_reads[1], "HIC2") == 0) return 0;
if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[1] == NULL)
{
fprintf(stderr, "[ERROR] lack r2 of HiC reads (--h2)\n");
return 0;
}
if(asm_opt->hic_files[1] != NULL && asm_opt->hic_files[0] == NULL)
if(asm_opt->hic_reads[1] != NULL && asm_opt->hic_reads[0] == NULL)
{
fprintf(stderr, "[ERROR] lack r1 of HiC reads (--h1)\n");
return 0;
}
if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[1] != NULL &&
asm_opt->hic_reads[0]->n != asm_opt->hic_reads[1]->n)
{
fprintf(stderr, "[ERROR] wrong r1 and r2 of HiC reads (--h1 && --h2)\n");
return 0;
}
if(asm_opt->hic_enzymes != NULL && asm_opt->hic_enzymes->n == 0)
{
fprintf(stderr, "[ERROR] wrong HiC enzymes (--enzyme)\n");
return 0;
}
if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[0]->n == 0)
{
fprintf(stderr, "[ERROR] wrong r1 of HiC reads (--h1)\n");
return 0;
}
if(asm_opt->hic_reads[1] != NULL && asm_opt->hic_reads[1]->n == 0)
{
fprintf(stderr, "[ERROR] wrong r2 of HiC reads (--h2)\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);
@@ -423,58 +468,60 @@ void get_queries(int argc, char *argv[], ketopt_t* opt, hifiasm_opt_t* asm_opt)
}
}
void get_hic_enzymes(char *argv, hifiasm_opt_t* asm_opt)
void get_hic_enzymes(char *argv, enzyme** x, int check_name)
{
int i, k, pre_i, len = strlen(argv);
asm_opt->hic_enzymes = (enzyme*)calloc(1, sizeof(enzyme));
(*x) = (enzyme*)calloc(1, sizeof(enzyme));
if(len == 0)
{
asm_opt->hic_enzymes->n = 0;
asm_opt->hic_enzymes->l = NULL;
asm_opt->hic_enzymes->a = NULL;
(*x)->n = 0; (*x)->l = NULL; (*x)->a = NULL;
return;
}
asm_opt->hic_enzymes->n = 1;
(*x)->n = 1;
for (i = pre_i = 0; i < len; i++)
{
if(argv[i] == ',')
{
asm_opt->hic_enzymes->n++;
(*x)->n++;
continue;
}
if(argv[i] != 'A' && argv[i] != 'C' && argv[i] != 'G' && argv[i] != 'T' &&
argv[i] != 'a' && argv[i] != 'c' && argv[i] != 'g' && argv[i] != 't' &&
argv[i] != 'N' && argv[i] != 'n')
if(check_name)
{
asm_opt->hic_enzymes->n = 0;
asm_opt->hic_enzymes->l = NULL;
asm_opt->hic_enzymes->a = NULL;
return;
if(argv[i] != 'A' && argv[i] != 'C' && argv[i] != 'G' && argv[i] != 'T' &&
argv[i] != 'a' && argv[i] != 'c' && argv[i] != 'g' && argv[i] != 't' &&
argv[i] != 'N' && argv[i] != 'n')
{
(*x)->n = 0;
(*x)->l = NULL;
(*x)->a = NULL;
return;
}
}
}
asm_opt->hic_enzymes->l = (int*)calloc(asm_opt->hic_enzymes->n, sizeof(int));
asm_opt->hic_enzymes->a = (char**)calloc(asm_opt->hic_enzymes->n, sizeof(char*));
(*x)->l = (int*)calloc((*x)->n, sizeof(int));
(*x)->a = (char**)calloc((*x)->n, sizeof(char*));
for (i = pre_i = k = 0; i < len; i++)
{
if(argv[i] == ',')
{
asm_opt->hic_enzymes->l[k] = i - pre_i;
asm_opt->hic_enzymes->a[k] = (char*)malloc(sizeof(char)*(asm_opt->hic_enzymes->l[k]+1));
memcpy(asm_opt->hic_enzymes->a[k], argv + pre_i, asm_opt->hic_enzymes->l[k]);
asm_opt->hic_enzymes->a[k][asm_opt->hic_enzymes->l[k]] = '\0';
(*x)->l[k] = i - pre_i;
(*x)->a[k] = (char*)malloc(sizeof(char)*((*x)->l[k]+1));
memcpy((*x)->a[k], argv + pre_i, (*x)->l[k]);
(*x)->a[k][(*x)->l[k]] = '\0';
pre_i = i + 1;
k++;
}
}
asm_opt->hic_enzymes->l[k] = i - pre_i;
asm_opt->hic_enzymes->a[k] = (char*)malloc(sizeof(char)*(asm_opt->hic_enzymes->l[k]+1));
memcpy(asm_opt->hic_enzymes->a[k], argv + pre_i, asm_opt->hic_enzymes->l[k]);
asm_opt->hic_enzymes->a[k][asm_opt->hic_enzymes->l[k]] = '\0';
(*x)->l[k] = i - pre_i;
(*x)->a[k] = (char*)malloc(sizeof(char)*((*x)->l[k]+1));
memcpy((*x)->a[k], argv + pre_i, (*x)->l[k]);
(*x)->a[k][(*x)->l[k]] = '\0';
}
int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
@@ -542,9 +589,10 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
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) asm_opt->hic_files[0] = opt.arg;
else if (c == 315) asm_opt->hic_files[1] = opt.arg;
else if (c == 316) get_hic_enzymes(opt.arg, asm_opt);
else if (c == 314) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[0]), 0);
else if (c == 315) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[1]), 0);
else if (c == 316) get_hic_enzymes(opt.arg, &(asm_opt->hic_enzymes), 1);
else if (c == 317) asm_opt->break_cov = 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);
@@ -563,7 +611,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
}
}
if (argc == opt.ind)
{
Print_H(asm_opt);
@@ -572,5 +620,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
get_queries(argc, argv, &opt, asm_opt);
return check_option(asm_opt);
}
CommandLines.h
+4 -3
View File
@@ -3,7 +3,7 @@
#include <pthread.h>
#define HA_VERSION "0.13-r308"
#define HA_VERSION "0.14-r309"
#define VERBOSE 0
@@ -35,7 +35,7 @@ typedef struct {
char *fn_bin_yak[2];
char *fn_bin_list[2];
char *extract_list;
char *hic_files[2];
enzyme *hic_reads[2];
enzyme *hic_enzymes;
int extract_iter;
int thread_num;
@@ -49,6 +49,7 @@ typedef struct {
double max_ov_diff_final;
int hom_cov;
int het_cov;
int break_cov;
int max_n_chain; // fall-back max number of chains to consider
int min_hist_kmer_cnt;
int load_index_from_disk;
@@ -106,7 +107,7 @@ static inline int ha_opt_triobin(const hifiasm_opt_t *opt)
static inline int ha_opt_hic(const hifiasm_opt_t *opt)
{
return ((opt->hic_files[0] && opt->hic_files[1]));
return ((opt->hic_reads[0] && opt->hic_reads[1]));
}
#endif
Overlaps.cpp
+1991 -493
View File
File diff suppressed because it is too large Load Diff
Overlaps.h
+10 -4
View File
@@ -21,6 +21,7 @@
#define PRIMARY_LABLE 0
#define ALTER_LABLE 1
#define HAP_LABLE 2
#define FAKE_LABLE 4
#define TRIO_THRES 0.9
#define DOUBLE_CHECK_THRES 0.1
#define FINAL_DOUBLE_CHECK_THRES 0.2
@@ -106,6 +107,11 @@ typedef struct {
uint8_t no_l_indel;
} asg_arc_t;
typedef struct {
size_t n, m;
asg_arc_t* a;
} kv_asg_arc_t;
typedef struct {
uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero
@@ -468,7 +474,7 @@ typedef struct {
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 set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8_t* flag);
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, int require_equal_nv, int test_tangle);
@@ -762,7 +768,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
{
qn = (node_max->a[b_max.readI]>>33);
set_R_to_U(ruIndex, qn, (b_max.b_0->b.a[b_max.untigI]>>1), 1);
set_R_to_U(ruIndex, qn, (b_max.b_0->b.a[b_max.untigI]>>1), 1, &(read_sg->seq[qn].c));
}
}
/*****************************label all unitigs****************************************/
@@ -819,7 +825,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
{
qn = (b_max.b_0->b.a[b_max.untigI]>>1);
set_R_to_U(ruIndex, qn, 1, 1);
set_R_to_U(ruIndex, qn, 1, 1, &(read_sg->seq[qn].c));
}
/*****************************label all reads****************************************/
@@ -1043,7 +1049,7 @@ typedef struct{
double weight;
uint32_t uID:31, del:1;
uint64_t dis;
uint64_t occ;
uint64_t occ:63, scaff:1;
///uint32_t enzyme;
} hc_edge;
Purge_Dups.cpp
+1 -1
View File
@@ -4278,7 +4278,7 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link)
for (i = 0, offset = 0; i < reads->n; i++)
{
rId = reads->a[i]>>33;
set_R_to_U(ruIndex, rId, uId, 1);
set_R_to_U(ruIndex, rId, uId, 1, &(read_g->seq[rId].c));
position_index[rId] = offset;
position_index[rId] = position_index[rId] << 32;
hic.cpp
+803 -378
View File
File diff suppressed because it is too large Load Diff
hic.h
+24 -1
View File
@@ -14,6 +14,21 @@ void push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t
void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link);
void hic_benchmark(ma_ug_t *ug, asg_t* read_g);
typedef struct {
double w;
uint32_t id, occ;
///uint32_t *bid, bid_n;
ma_utg_t *u;
uint64_t l_d, r_d;
}chain_hic_w_type;
typedef struct {
size_t n, m;
chain_hic_w_type* a;
uint32_t max_bub_id;
uint32_t *chain_idx, u_n;
}chain_hic_warp;
typedef struct {
long long g_occ, b_occ;
uint64_t id;
@@ -27,11 +42,12 @@ typedef struct {
kvec_t(uint32_t) num;
kvec_t(uint64_t) pathLen;
kvec_t(uint64_t) b_s_idx;
uint64_t s_bub, f_bub, b_bub, b_end_bub, tangle_bub, cross_bub;
uint64_t s_bub, f_bub, b_bub, b_end_bub, tangle_bub, cross_bub, mess_bub;
uint32_t check_het;
asg_t *b_g;
ma_ug_t* b_ug;
kvec_t(chain_w_type) chain_weight;
chain_hic_warp c_w;
} bubble_type;
#define P_het(B) ((B).num.n)
#define M_het(B) ((B).num.n + 1)
@@ -45,5 +61,12 @@ typedef struct {
void get_bubbles(bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t* sink, uint32_t** a, uint32_t* n, uint64_t* pathBase);
int load_hc_links(hc_links* link, const char *fn);
void write_hc_links(hc_links* link, const char *fn);
void destory_bubbles(bubble_type* bub);
void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link);
void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub);
uint32_t connect_bub_occ(bubble_type* bub, uint32_t root_id, uint32_t check_het);
void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, uint32_t check_het);
void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint32_t is_end);
void set_b_utg_weight_flag(bubble_type* bub, buf_t* b, uint32_t v, uint8_t* vis_flag, uint32_t flag, uint32_t* occ);
#endif