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

intergration correctio

Browse Source
This commit is contained in:
chhylp123
2022-03-01 20:25:53 -05:00
parent 0ce5a7ba5b
commit 9c6a3607d1
8 changed files with 1263 additions and 77 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.h
+1 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.16.2-r383"
#define HA_VERSION "0.16.2-r385"
#define VERBOSE 0
Correct.cpp
+48 -2
View File
@@ -8987,6 +8987,7 @@ inline void insert_snp_vv(haplotype_evdience_alloc* h, haplotype_evdience* a, ui
h->snp_stat.a[p->id].score = -1;
}
int insert_snp_ee(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, haplotype_evdience* u_a, UC_Read* g_read, void *km)
{
uint64_t i, m, occ_0, occ_1[5], occ_2, diff;
@@ -9003,6 +9004,7 @@ int insert_snp_ee(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a
// occ_2++;
// diff++;
// }
occ_2 += a[i].cov;
}
/**
@@ -9011,7 +9013,49 @@ int insert_snp_ee(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a
3. if occ_1 = 1, there are only one difference. It must be a sequencing error.
(for repeat, it maybe a snp at repeat. but ...)
**/
SnpStats *p = NULL;
uint32_t is_homopolymer = (uint32_t)-1;
if(occ_0 == 0 || diff <= 1) return 0;
for (i = m = 0; i < 4; i++) {
if(occ_1[i] >= 2){
if(!km) kv_pushp(SnpStats, h->snp_stat, &p);
else kv_pushp_km(km, SnpStats, h->snp_stat, &p);
p->id = h->snp_stat.n-1;
p->occ_0 = 1 + occ_0;
p->occ_1 = occ_1[i];
p->occ_2 = occ_2 - p->occ_0 - p->occ_1;
p->overlap_num = 0;
p->site = a[0].site;
p->score = -1;
p->overlap_num = occ_2;
if(is_homopolymer == (uint32_t)-1) {
is_homopolymer = if_is_homopolymer_strict(p->site, g_read->seq, g_read->length);
}
p->is_homopolymer = is_homopolymer;
occ_1[i] = p->id;
m++;
} else {
occ_1[i] = (uint64_t)-1;
}
}
if(m == 0) return 0;
for (i = m = 0; i < a_n; i++) {
if(a[i].type == 0) {
a[i].overlapSite = h->snp_stat.n-1;
} else if(occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]!=(uint64_t)-1){
a[i].cov = a[i].overlapSite;///note: only renew cov here!!!
a[i].overlapSite = occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]];
} else {
continue;
}
u_a[m++] = a[i];
}
/**
// if(c_snp && ovlp) {
// ;
// }
for (i = m = 0; i < 4; i++) {
if(occ_1[i] >= 2) {
insert_snp_vv(h, a, a_n, s_H[i], g_read, km);
@@ -9027,6 +9071,7 @@ int insert_snp_ee(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a
if(occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]] >= 2) u_a[m++] = a[i];
}
}
**/
return m;
}
@@ -9441,7 +9486,8 @@ 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,
int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec, void *km)
int force_repeat, int is_consensus, int* fully_cov, int* abnormal,
double max_ov_diff_ec, long long winLen, void *km)
{
clear_Correct_dumy(dumy, overlap_list, km);
@@ -9450,7 +9496,7 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref
Window_Pool w_inf;
init_Window_Pool(&w_inf, g_read->length, /**WINDOW_UL**//**WINDOW_UL_H**/MIN((((double)THRESHOLD_MAX_SIZE)/max_ov_diff_ec),WINDOW), (int)(1.0/max_ov_diff_ec));
init_Window_Pool(&w_inf, g_read->length, /**WINDOW_UL**//**WINDOW_UL_H**/winLen, (int)(1.0/max_ov_diff_ec));
int flag = 0;
while(get_Window(&w_inf, &window_start, &window_end) && flag != -2)
Correct.h
+2 -1
View File
@@ -1150,7 +1150,8 @@ 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,
int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec, void *km);
int force_repeat, int is_consensus, int* fully_cov, int* abnormal,
double max_ov_diff_ec, long long winLen, void *km);
/***
type:
Overlaps.cpp
+53 -19
View File
@@ -853,7 +853,7 @@ inline void set_reverse_overlap(ma_hit_t* dest, ma_hit_t* source)
dest->ml = source->ml;
dest->no_l_indel = source->no_l_indel;
/****************************may have bugs********************************/
dest->bl = Get_qe(*dest) - Get_qs(*dest);
dest->bl = source->bl/**Get_qe(*dest) - Get_qs(*dest)**/;
}
@@ -875,7 +875,7 @@ void normalize_ma_hit_t_single_side_advance(ma_hit_t_alloc* sources, long long n
qn = Get_qn(sources[i].buffer[j]);
tn = Get_tn(sources[i].buffer[j]);
sources[i].buffer[j].bl = Get_qe(sources[i].buffer[j]) - Get_qs(sources[i].buffer[j]);
// sources[i].buffer[j].bl = Get_qe(sources[i].buffer[j]) - Get_qs(sources[i].buffer[j]);
///if(sources[i].buffer[j].del) continue;
@@ -993,15 +993,15 @@ void normalize_ma_hit_t_single_side_aggressive(ma_hit_t_alloc* sources, long lon
qn = Get_qn(sources[i].buffer[j]);
tn = Get_tn(sources[i].buffer[j]);
sources[i].buffer[j].bl = Get_qe(sources[i].buffer[j]) - Get_qs(sources[i].buffer[j]);
// sources[i].buffer[j].bl = Get_qe(sources[i].buffer[j]) - Get_qs(sources[i].buffer[j]);
index = get_specific_overlap(&(sources[tn]), tn, qn);
if(index != -1)
{
sources[tn].buffer[index].bl = Get_qe(sources[tn].buffer[index])
- Get_qs(sources[tn].buffer[index]);
// sources[tn].buffer[index].bl = Get_qe(sources[tn].buffer[index])
// - Get_qs(sources[tn].buffer[index]);
if(Get_qs(sources[i].buffer[j]) == Get_ts(sources[tn].buffer[index])
&&
@@ -1124,15 +1124,15 @@ void normalize_ma_hit_t_single_side(ma_hit_t_alloc* sources, long long num_sourc
qn = Get_qn(sources[i].buffer[j]);
tn = Get_tn(sources[i].buffer[j]);
sources[i].buffer[j].bl = Get_qe(sources[i].buffer[j]) - Get_qs(sources[i].buffer[j]);
// sources[i].buffer[j].bl = Get_qe(sources[i].buffer[j]) - Get_qs(sources[i].buffer[j]);
index = get_specific_overlap(&(sources[tn]), tn, qn);
if(index != -1)
{
sources[tn].buffer[index].bl = Get_qe(sources[tn].buffer[index])
- Get_qs(sources[tn].buffer[index]);
// sources[tn].buffer[index].bl = Get_qe(sources[tn].buffer[index])
// - Get_qs(sources[tn].buffer[index]);
if(Get_qs(sources[i].buffer[j]) == Get_ts(sources[tn].buffer[index])
&&
@@ -9902,9 +9902,9 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
||
!check_weak_ma_hit_reverse(&(reverse_sources[qn]), sources, tn)**/)
{
sources[i].buffer[j].bl = 0;
sources[i].buffer[j].bl |= ((uint32_t)0x40000000);
index = get_specific_overlap(&(sources[tn]), tn, qn);
sources[tn].buffer[index].bl = 0;
sources[tn].buffer[index].bl |= ((uint32_t)0x40000000);
}
}
}
@@ -9919,13 +9919,14 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
{
if(sources[i].buffer[j].del) continue;
if(sources[i].buffer[j].bl != 0)
if(sources[i].buffer[j].bl&((uint32_t)0x40000000))
{
sources[i].buffer[j].del = 0;
sources[i].buffer[j].del = 1;
sources[i].buffer[j].bl -= ((uint32_t)0x40000000);
}
else
{
sources[i].buffer[j].del = 1;
sources[i].buffer[j].del = 0;
}
}
}
@@ -31100,16 +31101,41 @@ char *get_outfile_name(char* output_file_name)
return buf;
}
void create_ul_info(ma_hit_t_alloc* sources, int max_hang, int min_ovlp, long long gap_fuzz)
asg_t *build_init_sg(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int64_t n_read,
int64_t min_dp, uint64_t* readLen, int64_t mini_overlap_length, int64_t max_hang_length,
ma_sub_t *coverage_cut, R_to_U* ruIndex)
{
asg_t *sg = NULL;
clean_weak_ma_hit_t(sources, reverse_sources, n_read);
///ma_hit_sub is just use to init coverage_cut,
///it seems we do not need ma_hit_cut & ma_hit_flt
ma_hit_sub(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut);
detect_chimeric_reads(sources, n_read, readLen, coverage_cut, asm_opt.max_ov_diff_final * 2.0);
ma_hit_cut(sources, n_read, readLen, mini_overlap_length, &coverage_cut);
///print_binned_reads(sources, n_read, coverage_cut);
ma_hit_flt(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
///fix_binned_reads(sources, n_read, coverage_cut);
///just need to deal with trio here
ma_hit_contained_advance(sources, n_read, coverage_cut, ruIndex, max_hang_length, mini_overlap_length);
sg = ma_sg_gen(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
return sg;
}
void create_ul_info(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int64_t max_hang, int64_t min_ovlp, int64_t gap_fuzz,
int64_t min_dp, uint64_t* readLen, ma_sub_t *coverage_cut, R_to_U* ruIndex)
{
ug_opt_t opt; memset(&opt, 0, sizeof(opt));
opt.sources = sources;
opt.reverse_sources = reverse_sources;
opt.max_hang = max_hang;
opt.min_ovlp = min_ovlp;
opt.gap_fuzz = gap_fuzz;
opt.min_dp = min_dp;
opt.readLen = readLen;
opt.coverage_cut = coverage_cut;
opt.ruIndex = ruIndex;
ul_load(&opt);
exit(1);
}
void clean_graph(
@@ -31145,10 +31171,18 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
{
memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t));
}
if(asm_opt.ar) create_ul_info(sources, max_hang_length, mini_overlap_length, gap_fuzz);
if(asm_opt.ar) {
create_ul_info(sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz,
min_dp, readLen, coverage_cut, ruIndex);
exit(1);
} else {
// sg = build_init_sg(sources, reverse_sources, n_read, min_dp, readLen, mini_overlap_length, max_hang_length,
// coverage_cut, ruIndex);
clean_weak_ma_hit_t(sources, reverse_sources, n_read);
}
///print_binned_reads(sources, n_read, coverage_cut);
clean_weak_ma_hit_t(sources, reverse_sources, n_read);
///ma_hit_sub is just use to init coverage_cut,
///it seems we do not need ma_hit_cut & ma_hit_flt
ma_hit_sub(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut);
Overlaps.h
+18 -2
View File
@@ -212,6 +212,10 @@ typedef struct {
uint32_t utg:31, ori:1, start, len;
} utg_intv_t;
typedef struct {
uint32_t x, s, e;
} utg_ct_t;
typedef struct {
uint32_t *idx;
kvec_t(uint64_t) interval;
@@ -219,13 +223,21 @@ typedef struct {
typedef struct {
kvec_t(uint64_t) idx;
kvec_t(uint32_t) rids;
kvec_t(utg_ct_t) rids;
kvec_t(uint8_t) is_c;
} ul_contain;
typedef struct {
ma_ug_t *ug;
asg_t *rg;
uint64_t *idx;
} cvert_t;
typedef struct {
ma_ug_t *ug;
ucov_t *cc;
ul_contain *ct;
// cvert_t *nug;
// kv_ul_ov_t *ov;
} ul_idx_t;
@@ -985,7 +997,9 @@ typedef struct{
int min_ovlp;
int is_bench;
long long gap_fuzz;
int64_t min_dp;
bub_label_t* b_mask_t;
uint64_t* readLen;
}ug_opt_t;
void adjust_utg_by_trio(ma_ug_t **ug, asg_t* read_g, uint8_t flag, float drop_rate,
@@ -1034,7 +1048,9 @@ uint32_t tn, kv_u_trans_hit_t* ktb, uint32_t bn);
void clean_u_trans_t_idx(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g);
uint32_t test_dbug(ma_ug_t* ug, FILE* fp);
void write_dbug(ma_ug_t* ug, FILE* fp);
asg_t *build_init_sg(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int64_t n_read,
int64_t min_dp, uint64_t* readLen, int64_t mini_overlap_length, int64_t max_hang_length,
ma_sub_t *coverage_cut, R_to_U* ruIndex);
#define JUNK_COV 5
#define DISCARD_RATE 0.8
Process_Read.cpp
+132 -4
View File
@@ -7,6 +7,7 @@
#include "Correct.h"
#include "kalloc.h"
#define UL_FLANK 512
uint8_t seq_nt6_table[256] = {
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
@@ -749,7 +750,7 @@ void destory_Debug_reads(Debug_reads* x)
void init_all_ul_t(all_ul_t *x, All_reads *hR) {
memset(x, 0, sizeof(*x));
x->hR = hR; x->mm = 0x7fffffff;
x->hR = hR; x->mm = 0x40000000;
init_aux_table();
}
void destory_all_ul_t(all_ul_t *x) {
@@ -761,6 +762,7 @@ void destory_all_ul_t(all_ul_t *x) {
for (i = 0; i < x->nid.n; i++) free(x->nid.a[i].a);
free(x->nid.a);
free(x->ridx.idx.a); free(x->ridx.occ.a);
}
void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn, uint64_t nn_offset)
@@ -843,7 +845,7 @@ void push_subblock_original_bases(char* str, all_ul_t *x, ul_vec_t *p, uint32_t
}
}
void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on) {
void append_ul_t_compress_ovlp(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on) {
int64_t i, mine, maxs, ovlp, end;
ul_vec_t *p = NULL;
nid_t *np = NULL;
@@ -926,6 +928,132 @@ void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str,
}
}
void debug_append_ul_t(ul_ov_t *o, int64_t on, ul_vec_t *p)
{
int64_t k, l = 0;
uint32_t sp = (uint32_t)-1, ep = (uint32_t)-1, qss, qee, m;
for (k = on-1; k >= 0; k--) {
if(sp == (uint32_t)-1 || o[k].qe <= sp) {
if(sp != (uint32_t)-1) l += ep - sp;
if(sp == (uint32_t)-1) qss = o[k].qe, qee = p->rlen;
else qss = o[k].qe, qee = sp;
if(qee > qss) {
for (m = 0; m < p->bb.n; m++) {
if(qss == (p->bb.a[m].qs+((p->bb.a[m].hid>>15)&(0x7fffU))) &&
qee == (p->bb.a[m].qe-(p->bb.a[m].hid&(0x7fffU)))) {
break;
}
}
if(m >= p->bb.n) fprintf(stderr, "ERROR\n");
}
sp = o[k].qs;
ep = o[k].qe;
} else {
sp = MIN(sp, o[k].qs);
}
}
if(sp != (uint32_t)-1) l += ep - sp;
if(sp == (uint32_t)-1) qss = 0, qee = p->rlen;
else qss = 0, qee = sp;
if(qee > qss) {
for (m = 0; m < p->bb.n; m++) {
if(qss == (p->bb.a[m].qs+((p->bb.a[m].hid>>15)&(0x7fffU))) &&
qee == (p->bb.a[m].qe-(p->bb.a[m].hid&(0x7fffU)))) {
break;
}
}
if(m >= p->bb.n) fprintf(stderr, "ERROR\n");
}
}
void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on) {
int64_t i, mine, maxs, ovlp, st, et;
uint32_t o_l, o_r;
ul_vec_t *p = NULL;
nid_t *np = NULL;
ul_ov_t *z = NULL;
uc_block_t *b = NULL;
if(id) {
kv_pushp(nid_t, x->nid, &np);
np->n = id_l; MALLOC(np->a, np->n+1); memcpy(np->a, id, id_l); np->a[id_l] = '\0';
}
if(str) {
if(rid == NULL) {
kv_pushp(ul_vec_t, *x, &p);
memset(p, 0, sizeof(*p));
} else {
if((*rid) >= x->m) kv_resize(ul_vec_t, *x, (*rid) + 1);
if((*rid) >= x->n) {
memset(x->a+x->n, 0, sizeof(*p)*((*rid) + 1 - x->n));
x->n = (*rid) + 1;
}
p = &(x->a[(*rid)]);
}
p->bb.n = p->N_site.n = p->r_base.n = 0;
p->rlen = str_l;
if(o == NULL || on == 0) on = 0;
for (i = on-1, st = et = str_l; i >= 0; i--) {
z = &(o[i]);
mine = MIN(et, ((int64_t)z->qe)); maxs = MAX(st, ((int64_t)z->qs));
ovlp = mine - maxs;
if(ovlp < 0) {///push original bases
kv_pushp(uc_block_t, p->bb, &b);
b->hid = x->mm; b->rev = 0;
b->qe = maxs; b->qs = b->qe + ovlp;
o_l = (b->qs >= UL_FLANK?UL_FLANK:b->qs);
o_r = ((str_l-b->qe)>=UL_FLANK?UL_FLANK:(str_l-b->qe));
b->hid |= (o_l<<15); b->hid |= o_r;
b->qs -= o_l; b->qe += o_r;
b->ts = p->r_base.n; b->te = b->ts + B4L(b->qe-b->qs);
kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te;
ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site), b->qs);
}
///push ovlp bases
kv_pushp(uc_block_t, p->bb, &b);
b->hid = z->tn; b->rev = z->rev;
b->qs = z->qs; b->qe = z->qe;
b->ts = z->ts; b->te = z->te;
st = MIN(st, z->qs);
}
if(st > 0) {///push original bases
kv_pushp(uc_block_t, p->bb, &b);
b->hid = x->mm; b->rev = 0;
b->qe = st; b->qs = 0;
o_l = (b->qs >= UL_FLANK?UL_FLANK:b->qs);
o_r = ((str_l-b->qe)>=UL_FLANK?UL_FLANK:(str_l-b->qe));
b->hid |= (o_l<<15); b->hid |= o_r;
b->qs -= o_l; b->qe += o_r;
b->ts = p->r_base.n; b->te = b->ts + B4L(b->qe-b->qs);
kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te;
ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site), b->qs);
// push_subblock_original_bases(str, x, p, end, str_l, 321);//for debug
}
// debug_append_ul_t(o, on, p);
// char *sst = NULL; CALLOC(sst, str_l);//for debug
// retrieve_ul_t(NULL, sst, x, rid?*rid:x->n-1, 0, 0, -1);
// if(memcmp(sst, str, str_l)) {
// fprintf(stderr, "ap-Wrong read, id: %ld, [%d, %ld)\n", (int64_t)(rid?*rid:x->n-1), 0, str_l);
// for (i = 0; i < str_l; i++) {
// if(sst[i] != str[i]) fprintf(stderr, "[%ld] input:%c, decompress:%c\n", i, str[i], sst[i]);
// }
// }
// free(sst);
}
}
void retrieve_ul_t(UC_Read* i_r, char *i_s, all_ul_t *ref, uint64_t ID, uint8_t strand, int64_t s, int64_t l) {
ul_vec_t *p = &(ref->a[ID]);
if(l < 0) l = p->rlen;
@@ -956,7 +1084,7 @@ void retrieve_ul_t(UC_Read* i_r, char *i_s, all_ul_t *ref, uint64_t ID, uint8_t
sep = MIN(e, b->qe) - b->qs;
sl = sep - ssp;
if(b->hid == ref->mm){///original bases
if(b->hid&ref->mm){///original bases
offset = ssp&3;
begLen = 4-offset;
if(begLen > sl) begLen = sl;
@@ -1012,7 +1140,7 @@ void retrieve_ul_t(UC_Read* i_r, char *i_s, all_ul_t *ref, uint64_t ID, uint8_t
sl = sep - ssp;
///[ssp, sep)
if(b->hid == ref->mm){///original bases
if(b->hid&ref->mm){///original bases
begLen = sep&3;
offset = 4 - begLen;
if(begLen > sl) begLen = sl;
Process_Read.h
+6
View File
@@ -183,9 +183,15 @@ typedef struct
N_t N_site;
} ul_vec_t;
typedef struct{
kvec_t(uint32_t) idx;
kvec_t(uint64_t) occ;
} ul_vec_rid_t;
typedef struct
{
kvec_t(nid_t) nid;
ul_vec_rid_t ridx;
ul_vec_t *a;
size_t n, m;
All_reads *hR;
inter.cpp
+1003 -48
View File
File diff suppressed because it is too large Load Diff