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

variant calling

Browse Source
This commit is contained in:
chhylp123
2021-04-10 21:58:53 -04:00
parent 3218618bb4
commit 8f733750c4
5 changed files with 789 additions and 115 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Overlaps.cpp
+685 -98
View File
File diff suppressed because it is too large Load Diff
Overlaps.h
+13 -3
View File
@@ -52,6 +52,8 @@
#define CUT_DIF_HAP 12
///query is the read itself
typedef struct {
uint64_t qns;
@@ -113,15 +115,19 @@ typedef struct {
typedef struct {
uint32_t qs, qe, qn;
uint32_t ts, te, tn;
uint32_t nw;
uint8_t f:7, rev:1;
double nw;
uint8_t f:6, rev:1, del:1;
} u_trans_t;
typedef struct {
size_t n, m;
u_trans_t* a;
kvec_t(uint64_t) idx;
} kv_u_trans_t;
#define u_trans_a(x, id) ((x).a + ((x).idx.a[(id)]>>32))
#define u_trans_n(x, id) ((uint32_t)((x).idx.a[(id)]))
typedef struct {
uint64_t ul;
uint32_t v;
@@ -1144,7 +1150,7 @@ 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, uint32_t is_update_chain);
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, uint32_t is_update_chain);
uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes, hap_cov_t *cov, uint32_t is_update_chain, uint32_t keep_d);
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,
@@ -1176,6 +1182,10 @@ inline uint32_t get_origin_uid(uint32_t v, trans_chain* t_ch, uint32_t *off, uin
return (uint32_t)(((t_ch->rUidx[v>>1]>>1)<<1) + ((t_ch->rUidx[v>>1]^v)&1));
}
void get_chain_trans(trans_chain* t_ch, uint32_t id, uint32_t** x, uint32_t* x_occ, uint32_t** y, uint32_t* y_occ);
void chain_origin_trans_uid_by_distance(hap_cov_t *cov, asg_t *read_sg,
uint32_t *pri_a, uint32_t pri_n, uint32_t pri_beg, uint64_t *i_pri_len,
uint32_t *aux_a, uint32_t aux_n, uint32_t aux_beg, uint64_t *i_aux_len,
ma_ug_t *ug, uint32_t flag, double overall_score, const char* cmd);
#define JUNK_COV 5
#define DISCARD_RATE 0.8
Process_Read.h
+2
View File
@@ -102,6 +102,8 @@ typedef struct
#define NON_TRIO 4
#define DROP 5
#define SET_TRIO 8
#define CHAIN_MATCH 1
#define CHAIN_UNMATCH 0.334
typedef struct
{
Purge_Dups.cpp
+86 -11
View File
@@ -2331,7 +2331,7 @@ ma_hit_t_alloc* reverse_sources, long long xBegPos, long long xEndPos)
// fprintf(stderr, "tailIndex->a.n: %u, xReads->n: %u, total_match: %lld, hap_match: %lld, inp_match: %lld\n",
// tailIndex->a.n, xReads->n, (xEndPos - xBegPos + 1), hap_match, inp_match);
return ((double)(inp_match)*1) - ((double)(hap_match-inp_match)*0.334);
return ((double)(inp_match)*CHAIN_MATCH) - ((double)(hap_match-inp_match)*CHAIN_UNMATCH);
}
@@ -5368,6 +5368,88 @@ int max_hang, int min_ovlp, float drop_ratio, p_g_t *pg)
}
void chain_origin_trans_uid_by_purge(hap_overlaps *x, ma_ug_t *ug, hap_cov_t *cov, uint64_t* position_index)
{
uint32_t pri_uid, aux_uid, r_x, r_y;
hap_candidates hap_for, hap_rev, *hap = NULL;
long long x_pos_beg, x_pos_end, y_pos_beg, y_pos_end;
Get_rev(hap_for) = x->rev;
Get_x_beg(hap_for) = x->x_beg_id; Get_x_end(hap_for) = x->x_end_id - 1;
Get_y_beg(hap_for) = x->y_beg_id; Get_y_end(hap_for) = x->y_end_id - 1;
r_x = determine_hap_overlap_type_advance(&hap_for, &(ug->u.a[x->xUid]), &(ug->u.a[x->yUid]),
cov->ruIndex, cov->reverse_sources, cov->coverage_cut, cov->read_g, position_index,
cov->max_hang, cov->min_ovlp, x->xUid, x->yUid, &(cov->u_buffer), &(cov->tailIndex),
&(cov->prevIndex), &x_pos_beg, &x_pos_end, &y_pos_beg, &y_pos_end);
Get_rev(hap_rev) = x->rev;
Get_x_beg(hap_rev) = x->y_beg_id; Get_x_end(hap_rev) = x->y_end_id - 1;
Get_y_beg(hap_rev) = x->x_beg_id; Get_y_end(hap_rev) = x->x_end_id - 1;
r_y = determine_hap_overlap_type_advance(&hap_rev, &(ug->u.a[x->yUid]), &(ug->u.a[x->xUid]),
cov->ruIndex, cov->reverse_sources, cov->coverage_cut, cov->read_g, position_index,
cov->max_hang, cov->min_ovlp, x->yUid, x->xUid, &(cov->u_buffer), &(cov->tailIndex),
&(cov->prevIndex), &x_pos_beg, &x_pos_end, &y_pos_beg, &y_pos_end);
if(r_x == (uint32_t)-1 && r_y == (uint32_t)-1)
{
fprintf(stderr, "ERROR\n");
return;
}
Get_rev(hap_for) = x->rev;
Get_x_beg(hap_for) = x->x_beg_id; Get_x_end(hap_for) = x->x_end_id - 1;
Get_y_beg(hap_for) = x->y_beg_id; Get_y_end(hap_for) = x->y_end_id - 1;
Get_rev(hap_rev) = x->rev;
Get_x_beg(hap_rev) = x->y_beg_id; Get_x_end(hap_rev) = x->y_end_id - 1;
Get_y_beg(hap_rev) = x->x_beg_id; Get_y_end(hap_rev) = x->x_end_id - 1;
if(r_x != (uint32_t)-1 && r_y == (uint32_t)-1)
{
pri_uid = x->xUid; aux_uid = x->yUid; hap = &hap_for;
}
else if(r_x == (uint32_t)-1 && r_y != (uint32_t)-1)
{
aux_uid = x->xUid; pri_uid = x->yUid; hap = &hap_rev;
}
else
{
if(hap_for.score >= hap_rev.score)
{
pri_uid = x->xUid; aux_uid = x->yUid; hap = &hap_for;
}
else
{
aux_uid = x->xUid; pri_uid = x->yUid; hap = &hap_rev;
}
}
determine_hap_overlap_type_advance(hap, &(ug->u.a[pri_uid]), &(ug->u.a[aux_uid]),
cov->ruIndex, cov->reverse_sources, cov->coverage_cut, cov->read_g, position_index,
cov->max_hang, cov->min_ovlp, pri_uid, aux_uid, &(cov->u_buffer), &(cov->tailIndex),
&(cov->prevIndex), &x_pos_beg, &x_pos_end, &y_pos_beg, &y_pos_end);
uint64_t pri_len = ug->u.a[pri_uid].len, aux_len = ug->u.a[aux_uid].len;
pri_uid <<= 1; aux_uid <<= 1; aux_uid += hap->rev;
// uint32_t i_n = cov->t_ch->k_trans.n, i;
chain_origin_trans_uid_by_distance(cov, cov->read_g, &pri_uid, 1, x_pos_beg, &pri_len,
&aux_uid, 1, y_pos_beg, &aux_len, ug, RC_2, hap->score, __func__);
// fprintf(stderr, "\nocc: %u\n", (uint32_t)(cov->t_ch->k_trans.n - i_n));
// fprintf(stderr, "#s-utg%.6ul\t%u\t%u\td-utg%.6ul\t%u\t%u\trev(%u)\n",
// x->xUid+1, x->x_beg_pos, x->x_end_pos, x->yUid+1, x->y_beg_pos, x->y_end_pos, x->rev);
// for (i = i_n; i < cov->t_ch->k_trans.n; i++)
// {
// fprintf(stderr, "s-utg%.6ul\t%u\t%u\td-utg%.6ul\t%u\t%u\trev(%u)\n",
// cov->t_ch->k_trans.a[i].qn+1, cov->t_ch->k_trans.a[i].qs, cov->t_ch->k_trans.a[i].qe,
// cov->t_ch->k_trans.a[i].tn+1, cov->t_ch->k_trans.a[i].ts, cov->t_ch->k_trans.a[i].te,
// cov->t_ch->k_trans.a[i].rev);
// }
}
void collect_purge_trans_cov(ma_ug_t *ug, hap_overlaps_list* ha, hap_cov_t *cov, uint64_t* position_index)
{
uint32_t v, i, k, e, s, o, c_uId, p_uId, x_occ, y_occ;
@@ -5379,17 +5461,10 @@ void collect_purge_trans_cov(ma_ug_t *ug, hap_overlaps_list* ha, hap_cov_t *cov,
for (i = 0; i < ha->x[v].a.n; i++)
{
x = &(ha->x[v].a.a[i]);
/**
get_base_boundary_chain(cov->ruIndex, cov->reverse_sources, cov->coverage_cut,
cov->read_g, position_index, cov->max_hang, cov->min_ovlp, &(ug->u.a[x->xUid]),
&(ug->u.a[x->yUid]), x->xUid, x->yUid, x->x_beg_id, x->x_end_id-1,
x->y_beg_id, x->y_end_id-1, x->rev, &(cov->u_buffer), &(cov->tailIndex),
&(cov->prevIndex));
chain_origin_trans_uid(cov, cov->read_g, RC_2);
**/
if(x->yUid < x->xUid) continue;
chain_origin_trans_uid_by_purge(x, ug, cov, position_index);
q = &(ug->u.a[x->xUid]); s = x->x_beg_id; e = x->x_end_id; o = 0;
for (k = s, p_uId = (uint32_t)-1; k < e; k++)
{
hic.cpp
+3 -3
View File
@@ -2616,7 +2616,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag)
if(ug->g->seq[v>>1].del) continue;
if(asg_arc_n(ug->g, v) < 2) continue;
if((bub->index[v]&(uint32_t)3) != 0) continue;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0))
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0))
{
//beg is v, end is b.S.a[0]
//note b.b include end, does not include beg
@@ -2637,7 +2637,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag)
for (v = 0; v < n_vtx; ++v)
{
if((bub->index[v]&(uint32_t)3) !=2) continue;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0))
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0, 0))
{
//note b.b include end, does not include beg
i = b.b.n + 1;
@@ -2669,7 +2669,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag)
if((bub->num.a[k]>>31) == 0) bub->s_bub++;
v = (bub->num.a[k]<<1)>>1;
bub->num.a[k] = bub->list.n;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0))
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0, 0))
{
kv_push(uint64_t, bub->pathLen, pathLen);
//beg is v, end is b.S.a[0]