diff --git a/Overlaps.h b/Overlaps.h index 4a8a961..8a661cd 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -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;