mirror of
https://github.com/chhylp123/hifiasm.git
synced 2026-10-09 02:38:11 +08:00
for hifiasm 3.0
This commit is contained in:
chhylp123
+632
-21
@@ -1,5 +1,6 @@
|
||||
#ifndef __OVERLAPS__
|
||||
#define __OVERLAPS__
|
||||
#include <stdio.h>
|
||||
#include <stdint.h>
|
||||
#include "kvec.h"
|
||||
#include "kdq.h"
|
||||
@@ -20,6 +21,7 @@
|
||||
#define PRIMARY_LABLE 0
|
||||
#define ALTER_LABLE 1
|
||||
#define HAP_LABLE 2
|
||||
#define TRIO_THRES 0.9
|
||||
// #define PRIMARY_LABLE 1
|
||||
// #define ALTER_LABLE 2
|
||||
// #define HAP_LABLE 4
|
||||
@@ -40,6 +42,11 @@
|
||||
#define LONG_TIPS_UNDER_MAX_EXT 6
|
||||
#define LOOP 7
|
||||
|
||||
#define TRIM 10
|
||||
#define CUT 11
|
||||
#define CUT_DIF_HAP 12
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -65,6 +72,7 @@ typedef struct {
|
||||
void init_ma_hit_t_alloc(ma_hit_t_alloc* x);
|
||||
void clear_ma_hit_t_alloc(ma_hit_t_alloc* x);
|
||||
void resize_ma_hit_t_alloc(ma_hit_t_alloc* x, uint64_t size);
|
||||
void destory_all_ma_hit_t_alloc(ma_hit_t_alloc* x, uint64_t n_read);
|
||||
void destory_ma_hit_t_alloc(ma_hit_t_alloc* x);
|
||||
void add_ma_hit_t_alloc(ma_hit_t_alloc* x, ma_hit_t* element);
|
||||
void ma_hit_sort_tn(ma_hit_t *a, long long n);
|
||||
@@ -73,7 +81,6 @@ void ma_hit_sort_qns(ma_hit_t *a, long long n);
|
||||
int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources,
|
||||
char* output_file_name);
|
||||
|
||||
void normalize_ma_hit_t(ma_hit_t_alloc* sources, long long num_sources);
|
||||
|
||||
|
||||
typedef struct {
|
||||
@@ -99,6 +106,17 @@ typedef struct {
|
||||
uint8_t no_l_indel;
|
||||
} asg_arc_t;
|
||||
|
||||
|
||||
typedef struct {
|
||||
uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero
|
||||
uint32_t start, end; // start: starting vertex in the string graph; end: ending vertex
|
||||
uint32_t m, n; // number of reads
|
||||
uint64_t *a; // list of reads
|
||||
char *s; // unitig sequence is not null
|
||||
} ma_utg_t;
|
||||
|
||||
|
||||
|
||||
typedef struct {
|
||||
uint32_t len:31, del:1;
|
||||
uint8_t c;
|
||||
@@ -114,11 +132,26 @@ typedef struct {
|
||||
uint64_t *idx;
|
||||
|
||||
uint8_t* seq_vis;
|
||||
|
||||
uint32_t n_F_seq;
|
||||
ma_utg_t* F_seq;
|
||||
} asg_t;
|
||||
|
||||
typedef struct { size_t n, m; uint64_t *a; } asg64_v;
|
||||
|
||||
|
||||
typedef struct { size_t n, m; ma_utg_t *a; } ma_utg_v;
|
||||
|
||||
typedef struct {
|
||||
ma_utg_v u;
|
||||
asg_t *g;
|
||||
} ma_ug_t;
|
||||
|
||||
typedef struct {
|
||||
uint32_t utg:31, ori:1, start, len;
|
||||
} utg_intv_t;
|
||||
|
||||
|
||||
#define MA_HT_INT (-1)
|
||||
#define MA_HT_QCONT (-2)
|
||||
#define MA_HT_TCONT (-3)
|
||||
@@ -209,6 +242,23 @@ static inline int ma_hit2arc(const ma_hit_t *h, int ql, int tl, int max_hang, fl
|
||||
#define asg_arc_n(g, v) ((uint32_t)(g)->idx[(v)])
|
||||
#define asg_arc_a(g, v) (&(g)->arc[(g)->idx[(v)]>>32])
|
||||
|
||||
static inline uint32_t asg_get_arc(asg_t *g, uint32_t v, uint32_t w, asg_arc_t* t)
|
||||
{
|
||||
uint32_t i, nv = asg_arc_n(g, v);
|
||||
asg_arc_t *av = asg_arc_a(g, v);
|
||||
for (i = 0; i < nv; ++i)
|
||||
{
|
||||
if(av[i].del) continue;
|
||||
if(av[i].v == w)
|
||||
{
|
||||
(*t) = av[i];
|
||||
return 1;
|
||||
}
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
// append an arc
|
||||
static inline asg_arc_t *asg_arc_pushp(asg_t *g)
|
||||
{
|
||||
@@ -262,7 +312,8 @@ static inline void asg_seq_drop(asg_t *g, uint32_t s)
|
||||
///if output node is at primary
|
||||
/****************************may have hap bugs********************************/
|
||||
///if(g->seq[(av[i].v>>1)].c == PRIMARY_LABLE)
|
||||
if(g->seq[(av[i].v>>1)].c == PRIMARY_LABLE || g->seq[(av[i].v>>1)].c == HAP_LABLE)
|
||||
///if(g->seq[(av[i].v>>1)].c == PRIMARY_LABLE || g->seq[(av[i].v>>1)].c == HAP_LABLE)
|
||||
if(g->seq[(av[i].v>>1)].c != ALTER_LABLE)
|
||||
{/****************************may have hap bugs********************************/
|
||||
av[i].del = 1;
|
||||
asg_arc_del(g, av[i].v^1, v^1, 1);
|
||||
@@ -273,24 +324,7 @@ static inline void asg_seq_drop(asg_t *g, uint32_t s)
|
||||
}
|
||||
|
||||
|
||||
typedef struct {
|
||||
uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero
|
||||
uint32_t start, end; // start: starting vertex in the string graph; end: ending vertex
|
||||
uint32_t m, n; // number of reads
|
||||
uint64_t *a; // list of reads
|
||||
char *s; // unitig sequence is not null
|
||||
} ma_utg_t;
|
||||
|
||||
typedef struct { size_t n, m; ma_utg_t *a; } ma_utg_v;
|
||||
|
||||
typedef struct {
|
||||
ma_utg_v u;
|
||||
asg_t *g;
|
||||
} ma_ug_t;
|
||||
|
||||
typedef struct {
|
||||
uint32_t utg:31, ori:1, start, len;
|
||||
} utg_intv_t;
|
||||
|
||||
|
||||
/******************
|
||||
@@ -300,7 +334,10 @@ typedef struct {
|
||||
typedef struct {
|
||||
uint32_t p; // the optimal parent vertex
|
||||
uint32_t d; // the shortest distance from the initial vertex
|
||||
uint32_t c; // max count of reads
|
||||
uint32_t c; // max count of positive reads
|
||||
uint32_t m; // max count of negative reads
|
||||
uint32_t np; // max count of non-positive reads
|
||||
uint32_t nc; // max count of reads, no matter positive or negative
|
||||
uint32_t r:31, s:1; // r: the number of remaining incoming arc; s: state
|
||||
//s: state, s=0, this edge has not been visited, otherwise, s=1
|
||||
} binfo_t;
|
||||
@@ -331,6 +368,12 @@ typedef struct {
|
||||
uint64_t i;
|
||||
} kvec_t_u64_warp;
|
||||
|
||||
typedef struct {
|
||||
kvec_t(asg_arc_t) a;
|
||||
uint64_t i;
|
||||
}kvec_asg_arc_t_warp;
|
||||
|
||||
|
||||
|
||||
typedef struct {
|
||||
uint32_t q_pos;
|
||||
@@ -387,7 +430,6 @@ void add_overlaps(ma_hit_t_alloc* source_paf, ma_hit_t_alloc* dest_paf, uint64_t
|
||||
void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen);
|
||||
void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf,
|
||||
uint64_t* source_index, long long listLen);
|
||||
void print_revise_edges(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen);
|
||||
|
||||
#define EvaluateLen(U, id) ((U).a[(id)].start)
|
||||
#define IsMerge(U, id) ((U).a[(id)].end)
|
||||
@@ -421,5 +463,574 @@ void destory_R_to_U(R_to_U* x);
|
||||
void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig);
|
||||
void get_R_to_U(R_to_U* x, uint32_t rID, uint32_t* uID, uint32_t* is_Unitig);
|
||||
void transfor_R_to_U(R_to_U* x);
|
||||
void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, int require_equal_nv, int test_tangle);
|
||||
void clean_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
|
||||
long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold,
|
||||
R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen,
|
||||
uint32_t miniBiGraph, float chimeric_rate, int is_final_clean);
|
||||
int asg_pop_bubble_primary(asg_t *g, int max_dist);
|
||||
long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *g, long long circleLen, int is_drop);
|
||||
|
||||
typedef struct {
|
||||
asg_t* g;
|
||||
|
||||
asg_arc_t *av;
|
||||
uint32_t nv;
|
||||
uint32_t av_i;
|
||||
|
||||
asg_arc_t* new_edges;
|
||||
uint32_t new_edges_n;
|
||||
uint32_t new_edges_i;
|
||||
} Edge_iter;
|
||||
|
||||
void init_Edge_iter(asg_t* g, uint32_t v, asg_arc_t* new_edges, uint32_t new_edges_n, Edge_iter* x);
|
||||
int get_arc_t(Edge_iter* x, asg_arc_t* get);
|
||||
int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag);
|
||||
|
||||
|
||||
inline int get_real_length(asg_t *g, uint32_t v, uint32_t* v_s)
|
||||
{
|
||||
uint32_t i, kv = 0;
|
||||
for (i = 0, kv = 0; i < asg_arc_n(g, v); i++)
|
||||
{
|
||||
if(!asg_arc_a(g, v)[i].del)
|
||||
{
|
||||
if(v_s) v_s[kv] = asg_arc_a(g, v)[i].v;
|
||||
kv++;
|
||||
}
|
||||
}
|
||||
|
||||
return kv;
|
||||
}
|
||||
|
||||
inline uint32_t check_tip(asg_t *sg, uint32_t begNode, uint32_t* endNode, buf_t* b, uint32_t max_ext)
|
||||
{
|
||||
///cut tip of length <= max_ext
|
||||
uint32_t v = begNode, w;
|
||||
uint32_t kv;
|
||||
uint32_t eLen = 0;
|
||||
(*endNode) = (uint32_t)-1;
|
||||
b->b.n = 0;
|
||||
while (1)
|
||||
{
|
||||
kv = get_real_length(sg, v, NULL);
|
||||
(*endNode) = v;
|
||||
eLen++;
|
||||
if(b) kv_push(uint32_t, b->b, v);
|
||||
if(kv == 0) return END_TIPS;
|
||||
if(kv > 1) return MUL_OUTPUT;
|
||||
///if(eLen > max_ext) return LONG_TIPS;
|
||||
///kv must be 1 here
|
||||
kv = get_real_length(sg, v, &w);
|
||||
///here this value must be >= 1
|
||||
if(get_real_length(sg, w^1, NULL)!=1) return MUL_INPUT;
|
||||
v = w;
|
||||
if(v == begNode) return LOOP;
|
||||
if(eLen >= max_ext) return LONG_TIPS;
|
||||
}
|
||||
}
|
||||
|
||||
inline uint32_t get_unitig_back(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode,
|
||||
long long* nodeLen, long long* baseLen, buf_t* b)
|
||||
{
|
||||
ma_utg_v* u = NULL;
|
||||
uint32_t v = begNode, w, k;
|
||||
uint32_t kv;
|
||||
(*nodeLen) = (*baseLen) = 0;
|
||||
(*endNode) = (uint32_t)-1;
|
||||
if(ug!=NULL) u = &(ug->u);
|
||||
|
||||
while (1)
|
||||
{
|
||||
kv = get_real_length(sg, v, NULL);
|
||||
(*endNode) = v;
|
||||
if(u == NULL)
|
||||
{
|
||||
(*nodeLen)++;
|
||||
}
|
||||
else
|
||||
{
|
||||
(*nodeLen) += EvaluateLen((*u), v>>1);
|
||||
}
|
||||
if(b) kv_push(uint32_t, b->b, v);
|
||||
///means reach the end of a unitig
|
||||
if(kv!=1) (*baseLen) += sg->seq[v>>1].len;
|
||||
if(kv==0) return END_TIPS;
|
||||
if(kv>1) return MUL_OUTPUT;
|
||||
///kv must be 1 here
|
||||
kv = get_real_length(sg, v, &w);
|
||||
///means reach the end of a unitig
|
||||
if(get_real_length(sg, w^1, NULL)!=1)
|
||||
{
|
||||
(*baseLen) += sg->seq[v>>1].len;
|
||||
return MUL_INPUT;
|
||||
}
|
||||
|
||||
for (k = 0; k < asg_arc_n(sg, v); k++)
|
||||
{
|
||||
if(asg_arc_a(sg, v)[k].del) continue;
|
||||
///here is just one undeleted edge
|
||||
(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
|
||||
break;
|
||||
}
|
||||
|
||||
v = w;
|
||||
if(v == begNode) return LOOP;
|
||||
}
|
||||
}
|
||||
|
||||
inline uint32_t get_unitig(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode,
|
||||
long long* nodeLen, long long* baseLen, long long* max_stop_nodeLen, long long* max_stop_baseLen,
|
||||
uint32_t stops_threshold, buf_t* b)
|
||||
{
|
||||
ma_utg_v* u = NULL;
|
||||
uint32_t v = begNode, w, k;
|
||||
uint32_t kv, return_flag, n_stops = 0;
|
||||
long long pre_baseLen = 0, pre_nodeLen = 0;
|
||||
long long cur_baseLen = 0, cur_nodeLen = 0;
|
||||
(*max_stop_nodeLen) = (*max_stop_baseLen) = (*nodeLen) = (*baseLen) = 0;
|
||||
(*endNode) = (uint32_t)-1;
|
||||
if(ug!=NULL) u = &(ug->u);
|
||||
|
||||
while (1)
|
||||
{
|
||||
kv = get_real_length(sg, v, NULL);
|
||||
(*endNode) = v;
|
||||
if(u == NULL)
|
||||
{
|
||||
(*nodeLen)++;
|
||||
}
|
||||
else
|
||||
{
|
||||
(*nodeLen) += EvaluateLen((*u), v>>1);
|
||||
}
|
||||
if(b) kv_push(uint32_t, b->b, v);
|
||||
///means reach the end of a unitig
|
||||
if(kv!=1) (*baseLen) += sg->seq[v>>1].len;
|
||||
if(kv==0)
|
||||
{
|
||||
return_flag = END_TIPS;
|
||||
break;
|
||||
///return END_TIPS;
|
||||
}
|
||||
if(kv>1)
|
||||
{
|
||||
return_flag = MUL_OUTPUT;
|
||||
break;
|
||||
///return MUL_OUTPUT;
|
||||
}
|
||||
///kv must be 1 here
|
||||
kv = get_real_length(sg, v, &w);
|
||||
///means reach the end of a unitig
|
||||
if(get_real_length(sg, w^1, NULL)!=1)
|
||||
{
|
||||
|
||||
n_stops++;
|
||||
if(n_stops >= stops_threshold)
|
||||
{
|
||||
(*baseLen) += sg->seq[v>>1].len;
|
||||
return_flag = MUL_INPUT;
|
||||
break;
|
||||
///return MUL_INPUT;
|
||||
}
|
||||
else
|
||||
{
|
||||
for (k = 0; k < asg_arc_n(sg, v); k++)
|
||||
{
|
||||
if(asg_arc_a(sg, v)[k].del) continue;
|
||||
///here is just one undeleted edge
|
||||
(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
cur_baseLen = (*baseLen) - pre_baseLen;
|
||||
pre_baseLen = (*baseLen);
|
||||
if(cur_baseLen > (*max_stop_baseLen))
|
||||
{
|
||||
(*max_stop_baseLen) = cur_baseLen;
|
||||
}
|
||||
|
||||
|
||||
cur_nodeLen = (*nodeLen) - pre_nodeLen;
|
||||
pre_nodeLen = (*nodeLen);
|
||||
if(cur_nodeLen > (*max_stop_nodeLen))
|
||||
{
|
||||
(*max_stop_nodeLen) = cur_nodeLen;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
for (k = 0; k < asg_arc_n(sg, v); k++)
|
||||
{
|
||||
if(asg_arc_a(sg, v)[k].del) continue;
|
||||
///here is just one undeleted edge
|
||||
(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
v = w;
|
||||
if(v == begNode)
|
||||
{
|
||||
return_flag = LOOP;
|
||||
break;
|
||||
///return LOOP;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
cur_baseLen = (*baseLen) - pre_baseLen;
|
||||
pre_baseLen = (*baseLen);
|
||||
if(cur_baseLen > (*max_stop_baseLen))
|
||||
{
|
||||
(*max_stop_baseLen) = cur_baseLen;
|
||||
}
|
||||
|
||||
|
||||
cur_nodeLen = (*nodeLen) - pre_nodeLen;
|
||||
pre_nodeLen = (*nodeLen);
|
||||
if(cur_nodeLen > (*max_stop_nodeLen))
|
||||
{
|
||||
(*max_stop_nodeLen) = cur_nodeLen;
|
||||
}
|
||||
|
||||
return return_flag;
|
||||
}
|
||||
|
||||
#define UNAVAILABLE (uint32_t)-1
|
||||
#define PLOID 0
|
||||
#define NON_PLOID 1
|
||||
#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);
|
||||
}
|
||||
}
|
||||
/*****************************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;
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
/*****************************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);
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
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;
|
||||
uint32_t mother_occ;
|
||||
uint32_t ambig_occ;
|
||||
uint32_t drop_occ;
|
||||
uint32_t total;
|
||||
} Trio_counter;
|
||||
|
||||
void resolve_tangles(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig,
|
||||
long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t trio_flag,
|
||||
float drop_ratio);
|
||||
void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
|
||||
void rescue_contained_reads_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
|
||||
R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t chainLenThres, uint32_t is_bubble_check,
|
||||
uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes);
|
||||
void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
|
||||
R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t is_bubble_check,
|
||||
uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges);
|
||||
void deduplicate(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig,
|
||||
long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex);
|
||||
void all_to_all_deduplicate(ma_ug_t* ug, uint8_t postive_flag, float drop_rate,
|
||||
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
|
||||
void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
|
||||
void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
|
||||
ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges);
|
||||
void get_unitig_trio_flag(ma_utg_t* nsu, uint32_t flag, uint32_t* require, uint32_t* non_require, uint32_t* ambigious);
|
||||
void rescue_missing_overlaps_backward(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
|
||||
R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t backward_steps,
|
||||
uint32_t is_bubble_check, uint32_t is_primary_check);
|
||||
uint32_t get_edge_from_source(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
|
||||
R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t query, uint32_t target, asg_arc_t* t);
|
||||
|
||||
#endif
|
||||
Reference in New Issue
Block a user