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

smaller memory requirement

Browse Source
This commit is contained in:
chhylp123
2022-05-12 16:15:58 -04:00
parent dd8888aa90
commit a2198891e7
8 changed files with 1345 additions and 2581 deletions
Download Patch File Download Diff File Expand all files Collapse all files

43
Assembly.cpp
View File

@@ -461,6 +461,8 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug)
init_overlap_region_alloc(&b->olist);
init_overlap_region_alloc(&b->olist_hp);
init_fake_cigar(&(b->tmp_region.f_cigar));
memset(&(b->tmp_region.w_list), 0, sizeof(b->tmp_region.w_list));
CALLOC(b->tmp_region.w_list.a, 1); b->tmp_region.w_list.n = b->tmp_region.w_list.m = 1;
kv_init(b->b_buf.a);
kv_init(b->r_buf.a);
kv_init(b->k_flag.a);
@@ -488,6 +490,7 @@ void ha_ovec_destroy(ha_ovec_buf_t *b)
ha_abuf_destroy(b->ab);
ha_abufl_destroy(b->abl);
destory_fake_cigar(&(b->tmp_region.f_cigar));
free(b->tmp_region.w_list.a); free(b->tmp_region.w_list.c.a);
kv_destroy(b->b_buf.a);
kv_destroy(b->r_buf.a);
kv_destroy(b->k_flag.a);
@@ -529,21 +532,25 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b, int64_t *mem_a)
mem_clist += sizeof(*(b->k_flag.a.a)) * b->k_flag.a.m;
mem_clist += sizeof(*(b->sp.a)) * b->sp.m;
mem_clist += sizeof(*(b->tmp_region.f_cigar.buffer)) * b->tmp_region.f_cigar.size;
mem_clist += sizeof(*(b->tmp_region.w_list.a)) * b->tmp_region.w_list.n;
mem_clist += sizeof(*(b->tmp_region.w_list.c.a)) * b->tmp_region.w_list.c.n;
mem_olist = b->olist.size * sizeof(overlap_region);
for (i = 0; i < (int64_t)b->olist.size; ++i) {
const overlap_region *r = &b->olist.list[i];
mem_olist += r->w_list_size * sizeof(window_list);
mem_olist += (r->w_list.n*sizeof(*(r->w_list.a))) + (r->w_list.c.n*sizeof(*(r->w_list.c.a)));
mem_olist += r->f_cigar.size * 8;
mem_olist += r->boundary_cigars.size * sizeof(window_list);
mem_olist += (r->boundary_cigars.n*sizeof(*(r->boundary_cigars.a)))
+ (r->boundary_cigars.c.n*sizeof(*(r->boundary_cigars.c.a)));
}
mem_olist += b->olist_hp.size * sizeof(overlap_region);
for (i = 0; i < (int64_t)b->olist_hp.size; ++i) {
const overlap_region *r = &b->olist_hp.list[i];
mem_olist += r->w_list_size * sizeof(window_list);
mem_olist += (r->w_list.n*sizeof(*(r->w_list.a))) + (r->w_list.c.n*sizeof(*(r->w_list.c.a)));
mem_olist += r->f_cigar.size * 8;
mem_olist += r->boundary_cigars.size * sizeof(window_list);
mem_olist += (r->boundary_cigars.n*sizeof(*(r->boundary_cigars.a)))
+ (r->boundary_cigars.c.n*sizeof(*(r->boundary_cigars.c.a)));
}
if(b->ab) mem_ab += ha_abuf_mem(b->ab);
@@ -585,6 +592,13 @@ static void worker_ovec(void *data, long i, int tid)
{
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
int fully_cov, abnormal;
// if(i != 33) return;
// fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
// if (memcmp("m64012_190920_173625/88015004/ccs", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
// } else {
// return;
// }
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist,
0.02, asm_opt.max_n_chain, 1, NULL/**&(b->k_flag)**/, &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
@@ -593,7 +607,7 @@ static void worker_ovec(void *data, long i, int tid)
clear_Round2_alignment(&b->round2);
correct_overlap(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
&b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal);
&b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 1, &fully_cov, &abnormal);
b->num_read_base += b->self_read.length;
b->num_correct_base += b->correct.corrected_base;
@@ -624,6 +638,7 @@ static void worker_ovec(void *data, long i, int tid)
}
if(het_cnt) het_cnt[i] = get_het_cnt(&b->hap);
// fprintf(stderr, "[M::%s-end] rid->%ld\n", __func__, i);
}
@@ -653,7 +668,7 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
clear_Round2_alignment(&b->round2);
correct_overlap(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
&b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal);
&b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 1, &fully_cov, &abnormal);
b->num_read_base += b->self_read.length;
b->num_correct_base += b->correct.corrected_base;
@@ -919,10 +934,12 @@ void ha_overlap_and_correct(int round)
ha_opt_update_cov(&asm_opt, hom_cov);
het_cnt = NULL;
if(round == asm_opt.number_of_round-1 && asm_opt.is_dbg_het_cnt) CALLOC(het_cnt, R_INF.total_reads);
// fprintf(stderr, "[M::%s-start]\n", __func__);
if (asm_opt.required_read_name)
kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads);
else
kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads);///debug_for_fix
// fprintf(stderr, "[M::%s-end]\n", __func__);
if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
ha_pt_destroy(ha_idx);
@@ -1420,12 +1437,13 @@ void debug_affine_gap_alignment(overlap_region_alloc *overlap_list, UC_Read* g_r
kv_destroy(y_num);
}
/**
static void worker_ov_final_high_het(void *data, long i, int tid)
{
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, HIGH_HET_ERROR_RATE,
asm_opt.max_n_chain, 1, NULL/**&(b->k_flag)**/, &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
asm_opt.max_n_chain, 1, NULL, &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
overlap_region_sort_y_id(b->olist.list, b->olist.length);
ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length);
@@ -1449,6 +1467,7 @@ static void worker_ov_final_high_het(void *data, long i, int tid)
correct_overlap_high_het(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read);
push_final_overlaps_increment(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
}
**/
void Output_PAF()
{
@@ -1680,11 +1699,11 @@ void ha_overlap_final(void)
for (i = 0; i < asm_opt.thread_num; ++i)
b[i] = ha_ovec_init(asm_opt.flag & HA_F_HIGH_HET, 1,0);///b[i] = ha_ovec_init(1, 1);
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
if(asm_opt.flag & HA_F_HIGH_HET)
{
kt_for(asm_opt.thread_num, worker_ov_final_high_het, b, R_INF.total_reads);
}
else
// if(asm_opt.flag & HA_F_HIGH_HET)
// {
// kt_for(asm_opt.thread_num, worker_ov_final_high_het, b, R_INF.total_reads);
// }
// else
{
kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads);
}

3599
Correct.cpp
View File

File diff suppressed because it is too large Load Diff

42
Correct.h
View File

@@ -1017,28 +1017,6 @@ inline void destoryHaplotypeEvdience(haplotype_evdience_alloc* h)
destory_SNP_IDs(&(h->dp.SNP_IDs));
}
inline void destoryHaplotypeEvdience_buf(void *km, haplotype_evdience_alloc* h, int is_z)
{
kfree(km, h->list);
kfree(km, h->snp_stat.a);
kfree(km, h->snp_srt.a);
kfree(km, h->snp_matrix);
kfree(km, h->r_snp);
kfree(km, h->dp.backtrack);
kfree(km, h->dp.max);
kfree(km, h->dp.max_for_sort);
kfree(km, h->dp.visit);
kfree(km, h->dp.backtrack_length);
kfree(km, h->dp.buffer);
kfree(km, h->dp.max_buffer);
kfree(km, h->dp.SNP_IDs.buffer);
kfree(km, h->dp.SNP_IDs.IDs);
if(is_z) {
memset(h, 0, sizeof(*h));
memset(&(h->dp.SNP_IDs), 0, sizeof(h->dp.SNP_IDs));
h->dp.SNP_IDs.max_snp_id = -1;
}
}
inline void ResizeInitHaplotypeEvdience(haplotype_evdience_alloc* h)
{
@@ -1126,12 +1104,11 @@ void clear_Round2_alignment(Round2_alignment* h);
void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g, Graph* DAGCon,
Cigar_record* current_cigar, haplotype_evdience_alloc* hap,
Round2_alignment* second_round, int force_repeat, int is_consensus,
int* fully_cov, int* abnormal);
Round2_alignment* second_round, kvec_t_u64_warp* v_idx, window_list_alloc* win_ciagr_buf,
int force_repeat, int is_consensus, int* fully_cov, int* abnormal);
void init_Correct_dumy_buf(Correct_dumy* list, void *km);
void init_Correct_dumy(Correct_dumy* list);
void destory_Correct_dumy(Correct_dumy* list);
void destory_Correct_dumy_buf(void *km, Correct_dumy* list, int is_z);
void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list, void *km);
void clear_Correct_dumy_pure(Correct_dumy* list);
void get_seq_from_Graph(Graph* backbone, Graph* DAGCon, Correct_dumy* dumy, Cigar_record* current_cigar, char* self_string,
@@ -1150,6 +1127,7 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref
UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read,
Graph* g, Graph* DAGCon, Cigar_record* current_cigar,
haplotype_evdience_alloc* hap, Round2_alignment* second_round,
kvec_t_u64_warp* v_idx, window_list_alloc* win_ciagr_buf,
int force_repeat, int is_consensus, int* fully_cov, int* abnormal,
double max_ov_diff_ec, long long winLen, void *km);
@@ -1253,8 +1231,8 @@ long long* max_q_pos, long long* max_t_pos, long long* global_score,
long long* extention_score, long long* q_boundary_score, long long* q_boundary_t_coordinate,
long long* t_boundary_score, long long* t_boundary_q_coordinate,
long long* droped, int mode);
void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_INF,
UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read);
// void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_INF,
// UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read);
long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* overlap_read, uint8_t* x_num,
uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen);
int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u8_warp* k_flag, float hp_rate, int rlen, FILE* fp);
@@ -1271,7 +1249,17 @@ inline int if_exact_match(char* x, long long xLen, char* y, long long yLen, long
return 0;
}
inline void get_cigar_cell(window_list *idx, window_list_alloc *cc, uint32_t i, uint8_t *c, uint32_t *len)
{
uint32_t p = cc->c.a[idx->cidx+i];
(*c) = (uint8_t)(p>>30); (*len) = (p<<2)>>2;
}
inline void push_cigar_cell(window_list_alloc *res, uint8_t c, uint32_t len)
{
uint32_t p = c; p <<= 30; p += len;
kv_push(uint32_t, res->c, p);
}
#define FORWARD_KSW 0
#define BACKWARD_KSW 1

91
Hash_Table.cpp
View File

@@ -24,9 +24,9 @@ void init_overlap_region_alloc(overlap_region_alloc* list)
///list->list = (overlap_region*)malloc(sizeof(overlap_region)*list->size);
list->list = (overlap_region*)calloc(list->size, sizeof(overlap_region));
uint64_t i;
for (i = 0; i < list->size; i++)
{
for (i = 0; i < list->size; i++) {
init_fake_cigar(&(list->list[i].f_cigar));
init_window_list_alloc(&(list->list[i].w_list));
init_window_list_alloc(&(list->list[i].boundary_cigars));
}
}
@@ -36,10 +36,9 @@ void clear_overlap_region_alloc(overlap_region_alloc* list)
list->length = 0;
list->mapped_overlaps_length = 0;
uint64_t i = 0;
for (i = 0; i < list->size; i++)
{
list->list[i].w_list_length = 0;
for (i = 0; i < list->size; i++) {
clear_fake_cigar(&(list->list[i].f_cigar));
clear_window_list_alloc(&(list->list[i].w_list));
clear_window_list_alloc(&(list->list[i].boundary_cigars));
}
}
@@ -47,30 +46,14 @@ void clear_overlap_region_alloc(overlap_region_alloc* list)
void destory_overlap_region_alloc(overlap_region_alloc* list)
{
uint64_t i = 0;
for (i = 0; i < list->size; i++)
{
if (list->list[i].w_list_size != 0)
{
free(list->list[i].w_list);
}
for (i = 0; i < list->size; i++) {
destory_fake_cigar(&(list->list[i].f_cigar));
destory_window_list_alloc(&(list->list[i].w_list));
destory_window_list_alloc(&(list->list[i].boundary_cigars));
}
free(list->list);
}
void destory_overlap_region_alloc_buf(void *km, overlap_region_alloc* list, int is_z)
{
uint64_t i = 0;
for (i = 0; i < list->size; i++) {
if(list->list[i].w_list_size>0) kfree(km, list->list[i].w_list);
if(list->list[i].f_cigar.size>0) kfree(km, list->list[i].f_cigar.buffer);
if(list->list[i].boundary_cigars.size>0) kfree(km, list->list[i].boundary_cigars.buffer);
}
kfree(km, list->list);
if(is_z) memset(list, 0, sizeof(*list));
}
int get_fake_gap_pos(Fake_Cigar* x, int index)
{
return (x->buffer[index]>>32);
@@ -910,32 +893,18 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r
void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error,
int extra_begin, int extra_end, int error_threshold, int blockLen, void *km)
{
long long length = region->x_pos_e - region->x_pos_s + 1;
///the length of window may large or small than WINDOW
/****************************may have bugs********************************/
uint64_t num_windows = length / blockLen + 4;
/****************************may have bugs********************************/
window_list *p = NULL;
kv_pushp(window_list, region->w_list, &p);
///w_list_length has alredy set to be 0 at clear_overlap_region_alloc
if (num_windows > region->w_list_size)
{
region->w_list_size = num_windows;
if(!km) REALLOC(region->w_list, region->w_list_size);
else KREALLOC(km, region->w_list, region->w_list_size);
}
region->w_list[region->w_list_length].x_start = x_start;
region->w_list[region->w_list_length].x_end = x_end;
region->w_list[region->w_list_length].y_start = y_start;
region->w_list[region->w_list_length].y_end = y_end;
region->w_list[region->w_list_length].error = error;
region->w_list[region->w_list_length].cigar.length = -1;
region->w_list[region->w_list_length].extra_begin = extra_begin;
region->w_list[region->w_list_length].extra_end = extra_end;
region->w_list[region->w_list_length].error_threshold = error_threshold;
region->w_list_length++;
p->x_start = x_start;
p->x_end = x_end;
p->y_start = y_start;
p->y_end = y_end;
p->error = error;
p->extra_begin = extra_begin;
p->extra_end = extra_end;
p->error_threshold = error_threshold;
p->cidx = p->clen = 0;
}
void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, uint64_t strand)
@@ -1118,33 +1087,23 @@ void resize_fake_cigar(Fake_Cigar* x, uint64_t size, void *km)
void init_window_list_alloc(window_list_alloc* x)
{
x->buffer = NULL;
x->length = 0;
x->size = 0;
memset(x, 0, sizeof((*x)));
}
void clear_window_list_alloc(window_list_alloc* x)
{
x->length = 0;
x->n = x->c.n = 0;
}
void destory_window_list_alloc(window_list_alloc* x)
{
if(x->size != 0)
{
free((x->buffer));
}
free(x->a); free(x->c.a);
}
void resize_window_list_alloc(window_list_alloc* x, long long size, void *km)
void resize_window_list_alloc(window_list_alloc* x, uint64_t size)
{
if(size > x->size){
x->size = size;
if(!km) REALLOC(x->buffer, x->size);
else KREALLOC(km, x->buffer, x->size);
}
long long i;
for (i = 0; i < x->size; i++) x->buffer[i].error = -1;
x->length = 0;
kv_resize(window_list, *x, size); x->n = x->c.n = 0;
uint64_t k;
for (k = 0; k < x->m; k++) x->a[k].error = -1;
x->c.n = 0;
}

43
Hash_Table.h
View File

@@ -44,32 +44,29 @@ typedef struct
uint64_t end_pos;
} k_mer_pos_list;
typedef struct
{
int C_L[CIGAR_MAX_LENGTH];
char C_C[CIGAR_MAX_LENGTH];
int length;
} CIGAR;
// typedef struct
// {
// int C_L[CIGAR_MAX_LENGTH];
// char C_C[CIGAR_MAX_LENGTH];
// int length;
// } CIGAR;
typedef struct
{
///the begining and end of a window, instead of the whole overlap
uint64_t x_start;
uint64_t x_end;
int y_end;
int y_start;
int extra_begin;
int extra_end;
int error_threshold;
int error;
CIGAR cigar;
int32_t x_start, x_end;
int32_t y_start, y_end;
int32_t extra_begin, extra_end;
int32_t error, error_threshold;
uint32_t cidx, clen;
//CIGAR cigar;
} window_list;
typedef struct
{
window_list* buffer;
int32_t length;
int32_t size;
size_t n, m;
window_list *a;
kvec_t(uint32_t) c;
} window_list_alloc;
typedef struct
@@ -100,11 +97,12 @@ typedef struct
int8_t strong;
uint32_t non_homopolymer_errors;
window_list* w_list;
uint32_t w_list_size;
uint32_t w_list_length;
// window_list* w_list;
// uint32_t w_list_size;
// uint32_t w_list_length;
Fake_Cigar f_cigar;
window_list_alloc w_list;
window_list_alloc boundary_cigars;
} overlap_region;
@@ -149,7 +147,6 @@ void destory_Candidates_list_buf(void *km, Candidates_list* l, int is_z);
void init_overlap_region_alloc(overlap_region_alloc* list);
void clear_overlap_region_alloc(overlap_region_alloc* list);
void destory_overlap_region_alloc(overlap_region_alloc* list);
void destory_overlap_region_alloc_buf(void *km, overlap_region_alloc* list, int is_z);
void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error,
int extra_begin, int extra_end, int error_threshold, int blockLen, void *km);
@@ -196,7 +193,7 @@ void resize_Chain_Data(Chain_Data* x, long long size, void *km);
void init_window_list_alloc(window_list_alloc* x);
void clear_window_list_alloc(window_list_alloc* x);
void destory_window_list_alloc(window_list_alloc* x);
void resize_window_list_alloc(window_list_alloc* x, long long size, void *km);
void resize_window_list_alloc(window_list_alloc* x, uint64_t size);
long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, double band_width_threshold, int max_skip, int x_readLen, int y_readLen, void *km);
int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp,
ma_utg_v *ua, int add_beg_end, void *km);

63
POA.cpp
View File

@@ -307,66 +307,53 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length,
}
void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_string, long long x_length,
char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end)
char* y_string, long long y_length, window_list *cigar_idx, window_list_alloc *cigar_s, long long backbone_start, long long backbone_end)
{
int x_i, y_i, cigar_i;
x_i = 0;
y_i = 0;
cigar_i = 0;
int operation;
int operationLen;
int i;
int last_operation = -1;
int64_t x_i = 0, y_i = 0, c_i = 0, c_n = cigar_idx->clen;
uint32_t i, operLen; uint8_t oper; int8_t last_oper = -1;
// if(currentNodeID == 366 && x_length == 9 && y_length == 9) {
// fprintf(stderr, "[M::%s] currentNodeID::%lld, x_length::%lld, c_n::%ld, cidx::%u, cigar_s_n::%lld\n", __func__,
// currentNodeID, x_length, c_n, cigar_idx->cidx, (long long)cigar_s->c.n);
// }
///note that node 0 is the start node
///0 is match, 1 is mismatch, 2 is up, 3 is left
///2 mean y has more bases, while 3 means x has more bases
while (cigar_i < cigar->length)
{
operation = cigar->C_C[cigar_i];
operationLen = cigar->C_L[cigar_i];
for (c_i = 0; c_i < c_n; c_i++) {
get_cigar_cell(cigar_idx, cigar_s, c_i, &oper, &operLen);
///match/mismatch
if (operation == 0 || operation == 1)
{
for (i = 0; i < operationLen; i++)
{
// if(currentNodeID == 366 && x_length == 9 && y_length == 9) {
// fprintf(stderr, "[M::%s] c_i::%ld, oper::%u, operLen::%u, last_oper::%d\n", __func__, c_i, oper, operLen, last_oper);
// }
if (oper == 0 || oper == 1) { ///match/mismatch
for (i = 0; i < operLen; i++) {
///if the previous node is insertion, this node might be mismatch/match
add_mismatchEdge_weight(backbone, currentNodeID, y_string[y_i], last_operation);
x_i++;
y_i++;
currentNodeID++;
add_mismatchEdge_weight(backbone, currentNodeID, y_string[y_i], last_oper);
x_i++; y_i++; currentNodeID++;
}
}///insertion
else if (operation == 2)
{
} else if (oper == 2) { ///insertion
///the begin and end of cigar cannot be 2, so -1 is right here
///if (operationLen <= CORRECT_INDEL_LENGTH)
{
add_insertionEdge_weight(backbone, currentNodeID, y_string + y_i, operationLen);
add_insertionEdge_weight(backbone, currentNodeID, y_string + y_i, operLen);
backbone->g_nodes.list[currentNodeID].num_insertions++;
}
y_i += operationLen;
}
else if (operation == 3)
{
y_i += operLen;
} else if (oper == 3) {
///3 means x has more bases, that means backbone has more bases
///like a mismatch (-)
///if (operationLen <= CORRECT_INDEL_LENGTH)
{
add_deletionEdge_weight(backbone, currentNodeID, operationLen);
add_deletionEdge_weight(backbone, currentNodeID, operLen);
}
currentNodeID += operationLen;
x_i += operationLen;
currentNodeID += operLen;
x_i += operLen;
}
last_operation = operation;
cigar_i++;
last_oper = oper;
}
}

2
POA.h
View File

@@ -431,7 +431,7 @@ uint64_t* get_Topo_Sort_Order(Node_alloc* list, int need_sort);
void init_Graph(Graph* g);
void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, long long* startID, long long* endID);
void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_string, long long x_length,
char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end);
char* y_string, long long y_length, window_list *cigar_idx, window_list_alloc *cigar_s, long long backbone_start, long long backbone_end);
void destory_Graph(Graph* g);
void clear_Graph(Graph* g);
void Perform_POA(Graph* g, overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read);

43
inter.cpp
View File

@@ -2295,9 +2295,9 @@ uint32_t overlap_statistics(overlap_region_alloc* olist, ma_ug_t *ug, int64_t *t
**/
if(tt){
uint32_t z;
for (z = 0; z < olist->list[k].w_list_length; z++) {
if(olist->list[k].w_list[z].y_end != -1) {
if(tt) *tt += olist->list[k].w_list[z].x_end+1-olist->list[k].w_list[z].x_start;
for (z = 0; z < olist->list[k].w_list.n; z++) {
if(olist->list[k].w_list.a[z].y_end != -1) {
if(tt) *tt += olist->list[k].w_list.a[z].x_end+1-olist->list[k].w_list.a[z].x_start;
}
}
}
@@ -2619,38 +2619,38 @@ double es_win_err(overlap_region* o, int64_t winLen, int64_t s, int64_t e)
os = (o->x_pos_s/winLen)*winLen;
si = (s-os)/winLen; ei = (e-os)/winLen;
for (k = si+1; k <= ei-1; k++) {
tLen += o->w_list[k].x_end+1-o->w_list[k].x_start;
if(o->w_list[k].y_end != -1) {
tErr += o->w_list[k].error;
tLen += o->w_list.a[k].x_end+1-o->w_list.a[k].x_start;
if(o->w_list.a[k].y_end != -1) {
tErr += o->w_list.a[k].error;
} else {
tErr += o->w_list[k].x_end+1-o->w_list[k].x_start;
tErr += o->w_list.a[k].x_end+1-o->w_list.a[k].x_start;
}
}
k = si;
maxS = MAX(s, (int64_t)(o->w_list[k].x_start)); minE = MIN(e, (int64_t)(o->w_list[k].x_end)) + 1;
maxS = MAX(s, (int64_t)(o->w_list.a[k].x_start)); minE = MIN(e, (int64_t)(o->w_list.a[k].x_end)) + 1;
ov = minE > maxS? minE - maxS:0;
if(ov == 0) {
fprintf(stderr, "WARNNING-1, o->w_list_length->%u, o->x_id->%u, s->%ld, e->%ld, w_list_s->%lu, w_list_e->%lu, winLen->%ld, o->x_pos_s->%u, o->x_pos_e->%u, si->%ld, flag->%d\n",
o->w_list_length, o->x_id, s, e, o->w_list[k].x_start, o->w_list[k].x_end, winLen, o->x_pos_s, o->x_pos_e, si, o->w_list[k].y_end);
fprintf(stderr, "WARNNING-1, o->w_list.n->%u, o->x_id->%u, s->%ld, e->%ld, w_list_s->%d, w_list_e->%d, winLen->%ld, o->x_pos_s->%u, o->x_pos_e->%u, si->%ld, flag->%d\n",
(uint32_t)o->w_list.n, o->x_id, s, e, o->w_list.a[k].x_start, o->w_list.a[k].x_end, winLen, o->x_pos_s, o->x_pos_e, si, o->w_list.a[k].y_end);
}
tLen += ov/**o->w_list[k].x_end+1-o->w_list[k].x_start**/;
if(o->w_list[k].y_end != -1) {
tErr += (ov*o->w_list[k].error)/(o->w_list[k].x_end+1-o->w_list[k].x_start);
if(o->w_list.a[k].y_end != -1) {
tErr += (ov*o->w_list.a[k].error)/(o->w_list.a[k].x_end+1-o->w_list.a[k].x_start);
} else {
tErr += ov/**o->w_list[k].x_end+1-o->w_list[k].x_start**/;
}
k = ei;
maxS = MAX(s, (int64_t)(o->w_list[k].x_start)); minE = MIN(e, (int64_t)(o->w_list[k].x_end)) + 1;
maxS = MAX(s, (int64_t)(o->w_list.a[k].x_start)); minE = MIN(e, (int64_t)(o->w_list.a[k].x_end)) + 1;
ov = minE > maxS? minE - maxS:0;
if(ov == 0) {
fprintf(stderr, "WARNNING-2, o->w_list_length->%u, o->x_id->%u, s->%ld, e->%ld, w_list_s->%lu, w_list_e->%lu, winLen->%ld, o->x_pos_s->%u, o->x_pos_e->%u, ei->%ld, flag->%d\n",
o->w_list_length, o->x_id, s, e, o->w_list[k].x_start, o->w_list[k].x_end, winLen, o->x_pos_s, o->x_pos_e, ei, o->w_list[k].y_end);
fprintf(stderr, "WARNNING-2, o->w_list.n->%u, o->x_id->%u, s->%ld, e->%ld, w_list_s->%d, w_list_e->%d, winLen->%ld, o->x_pos_s->%u, o->x_pos_e->%u, ei->%ld, flag->%d\n",
(uint32_t)o->w_list.n, o->x_id, s, e, o->w_list.a[k].x_start, o->w_list.a[k].x_end, winLen, o->x_pos_s, o->x_pos_e, ei, o->w_list.a[k].y_end);
}
tLen += ov/**o->w_list[k].x_end+1-o->w_list[k].x_start**/;
if(o->w_list[k].y_end != -1) {
tErr += (ov*o->w_list[k].error)/(o->w_list[k].x_end+1-o->w_list[k].x_start);
if(o->w_list.a[k].y_end != -1) {
tErr += (ov*o->w_list.a[k].error)/(o->w_list.a[k].x_end+1-o->w_list.a[k].x_start);
} else {
tErr += ov/**o->w_list[k].x_end+1-o->w_list[k].x_start**/;
}
@@ -5135,7 +5135,7 @@ static void worker_for_ul_scall_alignment(void *data, long i, int tid) // callba
b->self_read.seq = s->seq[i]; b->self_read.length = s->len[i]; b->self_read.size = 0;
correct_ul_overlap(&b->olist, s->uu, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
&b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul, winLen, km);
&b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul, winLen, km);
// uint64_t k;
// for (k = 0; k < b->olist.length; k++) {
@@ -5171,11 +5171,6 @@ static void worker_for_ul_scall_alignment(void *data, long i, int tid) // callba
//
// if(l1 == 0 && l2 > 0) fprintf(stderr, "[M::%s::%lu::no_match]\n", UL_INF.nid.a[s->id+i].a, s->len[i]);
// fprintf(stderr, "[M::%s::%lu::] l1->%u; l2->%u\n", UL_INF.nid.a[s->id+i].a, s->len[i], l1, l2);
if(km) {
destory_overlap_region_alloc_buf(km, &b->olist, 1);
destory_Correct_dumy_buf(km, &b->correct, 1);
destoryHaplotypeEvdience_buf(km, &b->hap, 1);
}
}
@@ -5203,7 +5198,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call
b->self_read.seq = s->seq[i]; b->self_read.length = s->len[i]; b->self_read.size = 0;
correct_ul_overlap(&b->olist, s->uu, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
&b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul, winLen, NULL);
&b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul, winLen, NULL);
// uint64_t k;
// for (k = 0; k < b->olist.length; k++) {