mirror of
https://github.com/chhylp123/hifiasm.git
synced 2026-09-24 12:38:11 +08:00
backpup
This commit is contained in:
chhylp123
+1
-292
@@ -710,302 +710,11 @@ uint32_t stops_threshold, buf_t* b)
|
||||
#define UNAVAILABLE (uint32_t)-1
|
||||
#define PLOID 0
|
||||
#define NON_PLOID 1
|
||||
#define DIFF_HAP_RATE 0.75
|
||||
// #define DIFF_HAP_RATE 0.75
|
||||
#define TRIO_DROP_THRES 0.9
|
||||
#define TRIO_DROP_LENGTH_THRES 0.8
|
||||
#define MAX_STOP_RATE 0.6
|
||||
#define TANGLE_MISSED_THRES 0.6
|
||||
///if ug == NULL, nsg should be equal to read_sg
|
||||
inline uint32_t check_different_haps(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg,
|
||||
uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1,
|
||||
R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
|
||||
{
|
||||
uint32_t vEnd, qn, tn, j, is_Unitig, uId;
|
||||
long long ELen_0, ELen_1, tmp, max_stop_nodeLen, max_stop_baseLen;
|
||||
|
||||
b_0->b.n = b_1->b.n = 0;
|
||||
if(get_unitig(nsg, ug, v_0, &vEnd, &ELen_0, &tmp, &max_stop_nodeLen, &max_stop_baseLen,
|
||||
stops_threshold, b_0) == LOOP)
|
||||
{
|
||||
return UNAVAILABLE;
|
||||
}
|
||||
if(get_unitig(nsg, ug, v_1, &vEnd, &ELen_1, &tmp, &max_stop_nodeLen, &max_stop_baseLen,
|
||||
stops_threshold, b_1) == LOOP)
|
||||
{
|
||||
return UNAVAILABLE;
|
||||
}
|
||||
if(ELen_0<=min_edge_length || ELen_1<=min_edge_length) return UNAVAILABLE;
|
||||
|
||||
rIdContig b_max, b_min;
|
||||
b_max.b_0 = b_min.b_0 = NULL;
|
||||
b_max.offset = b_max.readI = b_max.untigI = 0;
|
||||
b_min.offset = b_min.readI = b_min.untigI = 0;
|
||||
|
||||
if(ELen_0<=ELen_1)
|
||||
{
|
||||
b_min.b_0 = b_0;
|
||||
b_max.b_0 = b_1;
|
||||
}
|
||||
else
|
||||
{
|
||||
b_min.b_0 = b_1;
|
||||
b_max.b_0 = b_0;
|
||||
}
|
||||
|
||||
uint32_t max_count = 0, min_count = 0;
|
||||
ma_utg_t *node_min = NULL, *node_max = NULL;
|
||||
if(ug != NULL)
|
||||
{
|
||||
/*****************************label all unitigs****************************************/
|
||||
for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
|
||||
{
|
||||
node_max = &(ug->u.a[b_max.b_0->b.a[b_max.untigI]>>1]);
|
||||
///each read
|
||||
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, &(read_sg->seq[qn].c));
|
||||
}
|
||||
}
|
||||
/*****************************label all unitigs****************************************/
|
||||
|
||||
///each unitig
|
||||
for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
|
||||
{
|
||||
|
||||
node_min = &(ug->u.a[(b_min.b_0->b.a[b_min.untigI]>>1)]);
|
||||
|
||||
///each read
|
||||
for (b_min.readI = 0; b_min.readI < node_min->n; b_min.readI++)
|
||||
{
|
||||
qn = node_min->a[b_min.readI]>>33;
|
||||
|
||||
/************************BUG: don't forget****************************/
|
||||
if(reverse_sources[qn].length > 0) min_count++;
|
||||
///if(reverse_sources[qn].length >= 0) min_count++;
|
||||
/************************BUG: don't forget****************************/
|
||||
for (j = 0; j < (long long)reverse_sources[qn].length; j++)
|
||||
{
|
||||
tn = Get_tn(reverse_sources[qn].buffer[j]);
|
||||
if(read_sg->seq[tn].del == 1)
|
||||
{
|
||||
get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
|
||||
if(tn == (uint32_t)-1 || is_Unitig == 1 || read_sg->seq[tn].del == 1) continue;
|
||||
}
|
||||
|
||||
get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
|
||||
if(uId!=(uint32_t)-1 && is_Unitig == 1)
|
||||
{
|
||||
max_count++;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
/*****************************label all unitigs****************************************/
|
||||
for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
|
||||
{
|
||||
node_max = &(ug->u.a[b_max.b_0->b.a[b_max.untigI]>>1]);
|
||||
///each read
|
||||
for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
|
||||
{
|
||||
qn = (node_max->a[b_max.readI]>>33);
|
||||
ruIndex->index[qn] = (uint32_t)-1;
|
||||
}
|
||||
}
|
||||
/*****************************label all unitigs****************************************/
|
||||
}
|
||||
else
|
||||
{
|
||||
/*****************************label all reads****************************************/
|
||||
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, &(read_sg->seq[qn].c));
|
||||
}
|
||||
/*****************************label all reads****************************************/
|
||||
|
||||
///each read
|
||||
for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
|
||||
{
|
||||
qn = (b_min.b_0->b.a[b_min.untigI]>>1);
|
||||
|
||||
/************************BUG: don't forget****************************/
|
||||
if(reverse_sources[qn].length > 0) min_count++;
|
||||
///if(reverse_sources[qn].length >= 0) min_count++;
|
||||
/************************BUG: don't forget****************************/
|
||||
|
||||
for (j = 0; j < (long long)reverse_sources[qn].length; j++)
|
||||
{
|
||||
tn = Get_tn(reverse_sources[qn].buffer[j]);
|
||||
if(nsg->seq[tn].del == 1)
|
||||
{
|
||||
get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
|
||||
if(tn == (uint32_t)-1 || is_Unitig == 1 || nsg->seq[tn].del == 1) continue;
|
||||
}
|
||||
|
||||
|
||||
get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
|
||||
if(uId!=(uint32_t)-1 && is_Unitig == 1)
|
||||
{
|
||||
max_count++;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*****************************label all reads****************************************/
|
||||
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);
|
||||
ruIndex->index[qn] = (uint32_t)-1;
|
||||
}
|
||||
/*****************************label all reads****************************************/
|
||||
}
|
||||
|
||||
// if(((v_0==7707) && (v_1==26867))||((v_1==7707) && (v_0==26867)))
|
||||
// {
|
||||
// fprintf(stderr, "******\nv_0>>1: %u, v_0&1: %u, ELen_0: %u\n", v_0>>1, v_0&1, (uint32_t)ELen_0);
|
||||
// fprintf(stderr, "v_1>>1: %u, v_1&1: %u, ELen_1: %u\n", v_1>>1, v_1&1, (uint32_t)ELen_1);
|
||||
// fprintf(stderr, "min_count: %u, max_count: %u, DIFF_HAP_RATE: %f\n\n",
|
||||
// min_count, max_count, DIFF_HAP_RATE);
|
||||
// }
|
||||
|
||||
if(min_count == 0) return UNAVAILABLE;
|
||||
if(max_count > min_count*DIFF_HAP_RATE) return PLOID;
|
||||
return NON_PLOID;
|
||||
}
|
||||
|
||||
inline uint32_t check_different_haps_naive(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg,
|
||||
uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1,
|
||||
R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
|
||||
{
|
||||
uint32_t vEnd, qn, tn, j, is_Unitig;
|
||||
long long ELen_0, ELen_1, tmp, max_stop_nodeLen, max_stop_baseLen;
|
||||
|
||||
b_0->b.n = b_1->b.n = 0;
|
||||
if(get_unitig(nsg, ug, v_0, &vEnd, &ELen_0, &tmp, &max_stop_nodeLen, &max_stop_baseLen,
|
||||
stops_threshold, b_0) == LOOP)
|
||||
{
|
||||
return UNAVAILABLE;
|
||||
}
|
||||
if(get_unitig(nsg, ug, v_1, &vEnd, &ELen_1, &tmp, &max_stop_nodeLen, &max_stop_baseLen,
|
||||
stops_threshold, b_1) == LOOP)
|
||||
{
|
||||
return UNAVAILABLE;
|
||||
}
|
||||
|
||||
if(ELen_0<=min_edge_length || ELen_1<=min_edge_length) return UNAVAILABLE;
|
||||
|
||||
rIdContig b_max, b_min;
|
||||
b_max.b_0 = b_min.b_0 = NULL;
|
||||
b_max.offset = b_max.readI = b_max.untigI = 0;
|
||||
b_min.offset = b_min.readI = b_min.untigI = 0;
|
||||
|
||||
if(ELen_0<=ELen_1)
|
||||
{
|
||||
b_min.b_0 = b_0;
|
||||
b_max.b_0 = b_1;
|
||||
}
|
||||
else
|
||||
{
|
||||
b_min.b_0 = b_1;
|
||||
b_max.b_0 = b_0;
|
||||
}
|
||||
|
||||
uint32_t max_count = 0, min_count = 0;
|
||||
ma_utg_t *node_min = NULL, *node_max = NULL;
|
||||
|
||||
if(ug != NULL)
|
||||
{
|
||||
///each unitig
|
||||
for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
|
||||
{
|
||||
|
||||
node_min = &(ug->u.a[(b_min.b_0->b.a[b_min.untigI]>>1)]);
|
||||
|
||||
///each read
|
||||
for (b_min.readI = 0; b_min.readI < node_min->n; b_min.readI++)
|
||||
{
|
||||
qn = node_min->a[b_min.readI]>>33;
|
||||
|
||||
if(reverse_sources[qn].length > 0) min_count++;
|
||||
for (j = 0; j < (long long)reverse_sources[qn].length; j++)
|
||||
{
|
||||
tn = Get_tn(reverse_sources[qn].buffer[j]);
|
||||
if(read_sg->seq[tn].del == 1)
|
||||
{
|
||||
get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
|
||||
if(tn == (uint32_t)-1 || is_Unitig == 1 || read_sg->seq[tn].del == 1) continue;
|
||||
}
|
||||
|
||||
///each unitig
|
||||
for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
|
||||
{
|
||||
node_max = &(ug->u.a[b_max.b_0->b.a[b_max.untigI]>>1]);
|
||||
///each read
|
||||
for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
|
||||
{
|
||||
if(tn == (node_max->a[b_max.readI]>>33))
|
||||
{
|
||||
max_count++;
|
||||
goto end_check_different_haps_ug;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
end_check_different_haps_ug:;
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
///each read
|
||||
for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
|
||||
{
|
||||
qn = (b_min.b_0->b.a[b_min.untigI]>>1);
|
||||
|
||||
if(reverse_sources[qn].length > 0) min_count++;
|
||||
|
||||
for (j = 0; j < (long long)reverse_sources[qn].length; j++)
|
||||
{
|
||||
tn = Get_tn(reverse_sources[qn].buffer[j]);
|
||||
if(nsg->seq[tn].del == 1)
|
||||
{
|
||||
get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
|
||||
if(tn == (uint32_t)-1 || is_Unitig == 1 || nsg->seq[tn].del == 1) continue;
|
||||
}
|
||||
|
||||
///each read
|
||||
for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
|
||||
{
|
||||
if((b_max.b_0->b.a[b_max.untigI]>>1) == tn)
|
||||
{
|
||||
max_count++;
|
||||
goto end_check_different_haps_non_ug;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
end_check_different_haps_non_ug:;
|
||||
}
|
||||
}
|
||||
|
||||
// if(((v_0==7707) && (v_1==26867))||((v_1==7707) && (v_0==26867)))
|
||||
// {
|
||||
// fprintf(stderr, "******\nv_0>>1: %u, v_0&1: %u, ELen_0: %u\n", v_0>>1, v_0&1, (uint32_t)ELen_0);
|
||||
// fprintf(stderr, "v_1>>1: %u, v_1&1: %u, ELen_1: %u\n", v_1>>1, v_1&1, (uint32_t)ELen_1);
|
||||
// fprintf(stderr, "min_count: %u, max_count: %u, DIFF_HAP_RATE: %f\n\n",
|
||||
// min_count, max_count, DIFF_HAP_RATE);
|
||||
// }
|
||||
|
||||
if(min_count == 0) return UNAVAILABLE;
|
||||
if(max_count > min_count*DIFF_HAP_RATE) return PLOID;
|
||||
return NON_PLOID;
|
||||
}
|
||||
|
||||
|
||||
|
||||
typedef struct {
|
||||
uint32_t father_occ;
|
||||
|
||||
Reference in New Issue
Block a user