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

using coverage

Browse Source
This commit is contained in:
chhylp123
2020-07-12 16:03:01 -04:00
parent a25c37ed89
commit c01b593d7c
7 changed files with 522 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+2 -2
View File
@@ -1203,8 +1203,8 @@ int ha_assemble(void)
ha_extract_print_list(&R_INF, asm_opt.extract_iter, asm_opt.extract_list);
exit(0);
}
if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt);
///if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt);
if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt), ovlp_loaded = 2;
///if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt), ovlp_loaded = 2;
if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2;
CommandLines.cpp
+3
View File
@@ -117,6 +117,9 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->purge_level_trio = 0;
asm_opt->purge_simi_rate = 0.75;
asm_opt->purge_overlap_len = 1;
asm_opt->recover_atg_cov_min = -1;
asm_opt->recover_atg_cov_max = -1;
asm_opt->hom_global_coverage = -1;
}
void destory_opt(hifiasm_opt_t* asm_opt)
CommandLines.h
+4
View File
@@ -58,6 +58,9 @@ typedef struct {
int purge_level_primary;
int purge_level_trio;
int purge_overlap_len;
int recover_atg_cov_min;
int recover_atg_cov_max;
int hom_global_coverage;
float max_hang_rate;
float min_drop_rate;
@@ -71,6 +74,7 @@ typedef struct {
long long num_recorrected_bases;
long long mem_buf;
long long coverage;
} hifiasm_opt_t;
extern hifiasm_opt_t asm_opt;
Correct.h
+1
View File
@@ -5,6 +5,7 @@
#include "Levenshtein_distance.h"
#include "POA.h"
#include "Process_Read.h"
#include "Correct.h"
//#define CORRECT_THRESHOLD 0.70
#define CORRECT_THRESHOLD 0.60
Overlaps.cpp
+214 -4
View File
@@ -21231,7 +21231,7 @@ void unroll_simple_case_advance(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* reve
while (n_reduce > 0)
{
n_reduce = 0;
///break nearly circle, forget why...
n_reduce += asg_arc_del_simple_circle_untig(NULL, NULL, nsg, 100, 0);
/**
for (v = 0; v < n_vtx; ++v)
@@ -22039,6 +22039,205 @@ R_to_U* ruIndex)
}
void append_utg(ma_ug_t* ptg, ma_ug_t* atg)
{
uint64_t num_nodes = 0;
asg_t* nsg = atg->g;
uint32_t v, n_vtx = nsg->n_seq;
ma_utg_t *p;
for (v = 0; v < n_vtx; ++v)
{
if(nsg->seq[v].del || atg->u.a[v].m == 0) continue;
num_nodes++;
}
if(num_nodes == 0) return;
ptg->u.n = ptg->u.n + num_nodes;
if(ptg->u.n > ptg->u.m)
{
ptg->u.m = ptg->u.n;
ptg->u.a = (ma_utg_t*)realloc(ptg->u.a, ptg->u.m*sizeof(ma_utg_t));
}
ptg->u.n = ptg->u.n - num_nodes;
for (v = 0; v < atg->g->n_seq; ++v)
{
if(atg->g->seq[v].del || atg->u.a[v].m == 0) continue;
p = &(ptg->u.a[ptg->u.n]);
p->len = atg->u.a[v].len;
p->circ = atg->u.a[v].circ;
p->start = atg->u.a[v].start;
p->end = atg->u.a[v].end;
p->m = atg->u.a[v].m; atg->u.a[v].m = 0;
p->n = atg->u.a[v].n; atg->u.a[v].n = 0;
p->a = atg->u.a[v].a; atg->u.a[v].a = 0;
p->s = atg->u.a[v].s; atg->u.a[v].s = 0;
asg_seq_set(ptg->g, ptg->u.n, p->len, 0);
ptg->u.n++;
}
asg_cleanup(ptg->g);
}
void print_utg_coverage(ma_ug_t *ug, ma_sub_t* coverage_cut, uint32_t v, ma_hit_t_alloc* sources)
{
asg_t* nsg = ug->g;
uint32_t rId, k, j;
ma_utg_t* u = NULL;
ma_hit_t *h;
if(nsg->seq[v].del) return;
u = &(ug->u.a[v]);
if(u->m == 0) return;
long long R_bases = 0, C_bases = 0;
long long U_R_bases = 0, U_C_bases = 0;
for (k = 0; k < u->n; k++)
{
rId = u->a[k]>>33;
C_bases = 0;
R_bases = coverage_cut[rId].e - coverage_cut[rId].s;
for (j = 0; j < (uint64_t)(sources[rId].length); j++)
{
h = &(sources[rId].buffer[j]);
if(h->el != 1) continue;
C_bases += Get_qe((*h)) - Get_qs((*h));
}
U_R_bases += R_bases;
U_C_bases += C_bases;
C_bases = C_bases/R_bases;
fprintf(stderr, "%.*s\t%lld\n", (int)Get_NAME_LENGTH(R_INF, rId), Get_NAME(R_INF, rId), C_bases);
}
fprintf(stderr, "v: %u, coverage: %lld\n\n", v, U_C_bases/U_R_bases);
}
void recover_utg_by_coverage(ma_ug_t **ptg, asg_t* read_g, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, R_to_U* ruIndex)
{
if(asm_opt.recover_atg_cov_min == -1) return;
if(asm_opt.recover_atg_cov_max == -1) return;
if(asm_opt.recover_atg_cov_min > asm_opt.recover_atg_cov_max) return;
ma_ug_t *atg = NULL;
atg = ma_ug_gen_primary(read_g, ALTER_LABLE);
asg_t* nsg = atg->g;
uint32_t v, n_vtx = nsg->n_seq, k, j, rId, available_reads = 0, keep_atg = 0, tn, is_Unitig;
ma_utg_t* u = NULL;
ma_hit_t *h;
long long R_bases = 0, C_bases = 0, C_bases_primary = 0, C_bases_alter = 0;
for (v = 0; v < n_vtx; ++v)
{
if(nsg->seq[v].del) continue;
u = &(atg->u.a[v]);
if(u->m == 0) continue;
available_reads = 0;
for (k = 0; k < u->n; k++)
{
rId = u->a[k]>>33;
C_bases = C_bases_primary = C_bases_alter = 0;
R_bases = coverage_cut[rId].e - coverage_cut[rId].s;
for (j = 0; j < (uint64_t)(sources[rId].length); j++)
{
h = &(sources[rId].buffer[j]);
if(h->el != 1) continue;
tn = Get_tn((*h));
if(read_g->seq[tn].del == 1)
{
///get the id of read that contains it
get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
}
if(read_g->seq[tn].del == 1) continue;
if(read_g->seq[tn].c == ALTER_LABLE)
{
C_bases_alter += Get_qe((*h)) - Get_qs((*h));
}
else
{
C_bases_primary += Get_qe((*h)) - Get_qs((*h));
}
}
C_bases = C_bases_primary + C_bases_alter;
if(C_bases_alter < C_bases * ALTER_COV_THRES) continue;
C_bases = C_bases/R_bases;
if(C_bases >= asm_opt.recover_atg_cov_min && C_bases <= asm_opt.recover_atg_cov_max)
{
available_reads++;
}
}
if(available_reads < (u->n * 0.8) || available_reads == 0)
{
asg_seq_del(nsg, v);
if(u->m!=0)
{
u->m = u->n = 0;
free(u->a);
u->a = NULL;
}
}
else
{
///print_utg_coverage(atg, coverage_cut, v, sources);
///fprintf(stderr, "rId: %u\n", rId);
///fprintf(stderr, "%.*s\t%lld\n", (int)Get_NAME_LENGTH(R_INF, rId), Get_NAME(R_INF, rId), C_bases);
keep_atg++;
}
}
if(keep_atg > 0)
{
asg_cleanup(nsg);
asg_symm(nsg);
append_utg(*ptg, atg);
n_vtx = read_g->n_seq;
for (v = 0; v < n_vtx; v++)
{
read_g->seq[v].c = ALTER_LABLE;
}
nsg = (*ptg)->g;
n_vtx = nsg->n_seq;
for (v = 0; v < n_vtx; ++v)
{
if(nsg->seq[v].del) continue;
u = &((*ptg)->u.a[v]);
if(u->m == 0) continue;
for (k = 0; k < u->n; k++)
{
rId = u->a[k]>>33;
read_g->seq[rId].c = nsg->seq[v].c;
}
}
n_vtx = read_g->n_seq;
for (v = 0; v < n_vtx; v++)
{
if(read_g->seq[v].c == ALTER_LABLE)
{
asg_seq_drop(read_g, v);
}
}
}
fprintf(stderr, "keep_atg: %u\n", keep_atg);
ma_ug_destroy(atg);
}
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,
long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold,
@@ -22049,7 +22248,8 @@ kvec_asg_arc_t_warp* new_rtg_edges)
uint32_t v, n_vtx = nsg->n_seq, k, rId, just_contain;
ma_utg_t* u = NULL;
///print_utg_coverage(*ug, coverage_cut, 440, sources);
///exit(0);
/**
kvec_t_u32_warp new_rtg_nodes;
kv_init(new_rtg_nodes.a);
@@ -22110,6 +22310,7 @@ kvec_asg_arc_t_warp* new_rtg_edges)
renew_utg(ug, read_g, new_rtg_edges);
}
}
n_vtx = read_g->n_seq;
@@ -22134,7 +22335,6 @@ kvec_asg_arc_t_warp* new_rtg_edges)
}
}
n_vtx = read_g->n_seq;
for (v = 0; v < n_vtx; v++)
{
@@ -22143,7 +22343,17 @@ kvec_asg_arc_t_warp* new_rtg_edges)
asg_seq_drop(read_g, v);
}
}
if(asm_opt.recover_atg_cov_min == -1 || asm_opt.recover_atg_cov_max == -1)
{
asm_opt.recover_atg_cov_max = asm_opt.hom_global_coverage/HOM_PEAK_RATE;
asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.8;
asm_opt.recover_atg_cov_max = asm_opt.recover_atg_cov_max * 1.2;
}
fprintf(stderr, "asm_opt.recover_atg_cov_min: %d\n", asm_opt.recover_atg_cov_min);
fprintf(stderr, "asm_opt.recover_atg_cov_max: %d\n", asm_opt.recover_atg_cov_max);
recover_utg_by_coverage(ug, read_g, coverage_cut, sources, ruIndex);
/**
kv_destroy(new_rtg_nodes.a);
Purge_Dups.cpp
+293 -31
View File
@@ -6,6 +6,9 @@
#include "Overlaps.h"
#include "Correct.h"
#include "kthread.h"
#include "kdq.h"
KDQ_INIT(uint64_t)
#define Cal_Off(OFF) ((long long)((uint32_t)((OFF)>>32)) - (long long)((uint32_t)((OFF))))
#define Get_xOff(OFF) ((long long)((uint32_t)((OFF)>>32)))
@@ -149,7 +152,8 @@ void print_peak(long long* cov_buf, long long cov_buf_length, long long max_i)
}
void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo_peak_cov, long long* hom_peak, long long* het_peak)
void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo_peak_cov,
long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* coverage_only)
{
long long i, start, err_i, max_i, max2_i, max3_i, topo_peak_i, max, max2, max3, topo_peak, min;
@@ -240,36 +244,71 @@ void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo
{
topo_peak_i = (*topo_peak_cov);
topo_peak = cov_buf[topo_peak_i];
if(topo_peak_i <= max_i * 1.2 && topo_peak_i >= max_i * 0.8 && topo_peak > max * 0.05)
{
(*het_peak) = max_i;
}
///fprintf(stderr, "topo_peak: %lld, topo_peak_i: %lld\n", topo_peak, topo_peak_i);
if (topo_peak <= max * 0.05) topo_peak_i = topo_peak = -1;
}
long long k_mer_het, k_mer_hom, coverage_het, coverage_hom, alter_peak;
k_mer_het = k_mer_hom = coverage_het = coverage_hom = alter_peak = -1;
///if we really want to use rev_sources for double checking, we should use peak instead of mean
///if we found a small peak at the right hand of the largest peak
alter_peak = topo_peak_i;
k_mer_het = asm_opt.het_cov;
k_mer_hom = asm_opt.hom_cov;
if(max3_i > 0)
{
(*het_peak) = max_i;
(*hom_peak) = max3_i;
coverage_het = max_i;
coverage_hom = max3_i;
}
else if((*het_peak) == -1)
else
{
(*het_peak) = max2_i;
(*hom_peak) = max_i;
coverage_het = max2_i;
coverage_hom = max_i;
}
if(k_mer_het != -1)
{
(*het_peak) = k_mer_het;
(*hom_peak) = k_mer_hom;
return;
}
else if(coverage_het != -1)
{
(*het_peak) = coverage_het;
(*hom_peak) = coverage_hom;
return;
}
else if(k_mer_hom > coverage_hom*1.5)
{
(*het_peak) = coverage_hom;
(*hom_peak) = k_mer_hom;
return;
}
else if(alter_peak != -1)
{
///if peak is het, coverage peak is more reliable
if(coverage_hom >= alter_peak*0.8 && coverage_hom <= alter_peak*1.2)
{
(*het_peak) = coverage_hom;
return;
}///if peak is homo, k-mer peak is more reliable
else if(k_mer_hom >= alter_peak*0.8*2 && k_mer_hom <= alter_peak*1.2*2)
{
(*hom_peak) = k_mer_hom;
return;
}
}
(*k_mer_only) = k_mer_hom;
(*coverage_only) = coverage_hom;
// fprintf(stderr, "max: %lld, max_i: %lld\n", max, max_i);
// fprintf(stderr, "max2: %lld, max2_i: %lld\n", max2, max2_i);
// fprintf(stderr, "max3: %lld, max3_i: %lld\n", max3, max3_i);
fprintf(stderr, "[M::%s] Heterozygous k-mer peak: %d\n", __func__, asm_opt.het_cov);
fprintf(stderr, "[M::%s] Homozygous k-mer peak: %d\n", __func__, asm_opt.hom_cov);
fprintf(stderr, "[M::%s] Heterozygous coverage peak: %lld\n", __func__, (*het_peak));
fprintf(stderr, "[M::%s] Homozygous coverage peak: %lld\n", __func__, (*hom_peak));
fprintf(stderr, "[M::%s] Alter coverage peak: %lld\n", __func__, topo_peak_i);
// fprintf(stderr, "[M::%s] Heterozygous k-mer peak: %d\n", __func__, asm_opt.het_cov);
// fprintf(stderr, "[M::%s] Homozygous k-mer peak: %d\n", __func__, asm_opt.hom_cov);
// fprintf(stderr, "[M::%s] Heterozygous coverage peak: %lld\n", __func__, (*het_peak));
// fprintf(stderr, "[M::%s] Homozygous coverage peak: %lld\n", __func__, (*hom_peak));
// fprintf(stderr, "[M::%s] Alter coverage peak: %lld\n", __func__, topo_peak_i);
}
@@ -278,9 +317,7 @@ void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo
long long get_alter_peak(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, uint64_t* position_index,
ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
{
#define ALTER_COV_THRES 0.9
#define REAL_ALTER_THRES 0.1
ma_utg_t* u = NULL;
asg_t* nsg = ug->g;
uint64_t v, j, k, qn, n_vtx = nsg->n_seq, primary_bases = 0, alter_bases = 0;
@@ -372,7 +409,9 @@ ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
return max_i;
}
long long get_read_coverage_thres(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, uint64_t* position_index, ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, uint64_t n_read, long long cov_buf_length)
long long get_read_coverage_thres(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, uint64_t* position_index,
ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, uint64_t n_read, long long cov_buf_length,
long long* k_mer_only, long long* coverage_only)
{
uint64_t i, j;
long long* cov_buf = NULL;
@@ -401,13 +440,14 @@ long long get_read_coverage_thres(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, u
cov_buf_length);
}
get_read_peak(cov_buf, cov_buf_length, alter_peak == -1? NULL: &alter_peak, &hom_peak, &het_peak);
get_read_peak(cov_buf, cov_buf_length, alter_peak == -1? NULL: &alter_peak, &hom_peak, &het_peak,
k_mer_only, coverage_only);
free(cov_buf);
if(hom_peak != -1) return hom_peak*1.25;
if(het_peak != -1) return het_peak*2.50;
return cov_buf_length;
if(hom_peak != -1) return hom_peak*HOM_PEAK_RATE;
if(het_peak != -1) return het_peak*HET_PEAK_RATE;
return -1;
}
@@ -1901,7 +1941,7 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon
sources, coverage_cut);
///fprintf(stderr, "ploid_coverage: %lu, cov_threshold: %lu\n", ploid_coverage, cov_threshold);
if(ploid_coverage >= cov_threshold) return NON_PLOID;
if(cov_threshold > 0 && ploid_coverage >= cov_threshold) return NON_PLOID;
return PLOID;
}
@@ -2783,7 +2823,8 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid)
kvec_t_i32_warp* prevIndex_vec = &(hap_buf->buf[tid].u_buffer_prevIndex);
kvec_t_i32_warp* begIndex_vec = &(hap_buf->buf[tid].u_buffer_beg);
kvec_t_u8_warp* flag_vec = &(hap_buf->buf[tid].u_buffer_flag);
long long cov_threshold = hap_buf->cov_threshold;
uint64_t cov_threshold = hap_buf->cov_threshold;
if(hap_buf->cov_threshold < 0) cov_threshold = (uint64_t)-1;
ma_utg_t *xReads = NULL, *yReads = NULL;
ma_hit_t_alloc *xR = NULL;
@@ -3927,6 +3968,207 @@ void print_all_purge_ovlp(ma_ug_t *ug, hap_overlaps_list* all_ovlp)
}
}
inline int get_available_cnt(asg_t *g, uint32_t v, uint8_t* del, asg_arc_t* v_s)
{
//v has direction
if(del && del[v>>1]) return 0;
uint32_t i, kv = 0;
asg_arc_t *av = asg_arc_a(g, v);
uint32_t nv = asg_arc_n(g, v);
for (i = 0, kv = 0; i < nv; i++)
{
if(!av[i].del)
{
if(del && del[av[i].v>>1]) continue;
if(v_s) v_s[kv] = av[i];
kv++;
}
}
return kv;
}
long long get_specific_contig_length(asg_t *g, uint8_t *del)
{
asg_cleanup(g);
uint32_t v, n_vtx = g->n_seq * 2, q_occ;
uint8_t *mark = NULL;
///is a queue
//kdq_t(uint64_t) *q;
///each node has two directions
//q = kdq_init(uint64_t);
mark = (uint8_t*)calloc(n_vtx, 1);
long long totalLen = 0;
for (v = 0; v < n_vtx; ++v)
{
uint32_t w, x, l, start, end, len;
asg_arc_t arc;
if (g->seq[v>>1].del || mark[v]) continue;
if (get_available_cnt(g, v, del, NULL) == 0 && get_available_cnt(g, (v^1), del, NULL) != 0) continue;
if (del[v>>1]) continue;
mark[v] = 1;
//q->count = 0, start = v, end = v^1, len = 0;
q_occ =0, start = v, end = v^1, len = 0;
// forward
w = v;
while (1)
{
/**
* w----->x
* w<-----x
* that means the only suffix of w is x, and the only prefix of x is w
**/
if (get_available_cnt(g, w, del, NULL) != 1) break;
get_available_cnt(g, w, del, &arc);
x = arc.v; // w->x
if (get_available_cnt(g, x^1, del, NULL) != 1) break;
/**
* another direction of w would be marked as used (since w has been used)
**/
mark[x] = mark[w^1] = 1;
///l is the edge length, instead of overlap length
///note: edge length is different with overlap length
///l = asg_arc_len(arc_first(g, w));
get_available_cnt(g, w, del, &arc);
l = ((uint32_t)((arc).ul));
//kdq_push(uint64_t, q, (uint64_t)w<<32 | l);
q_occ++;
end = x^1, len += l;
w = x;
if (x == v) break;
}
//if (start != (end^1) || kdq_size(q) == 0) { // linear unitig
if (start != (end^1) || q_occ == 0) { // linear unitig
///length of seq, instead of edge
l = g->seq[end>>1].len;
//kdq_push(uint64_t, q, (uint64_t)(end^1)<<32 | l);
q_occ++;
len += l;
} else { // circular unitig
start = end = UINT32_MAX;
goto add_unitig; // then it is not necessary to do the backward
}
// backward
x = v;
while (1) { // similar to forward but not the same
if (get_available_cnt(g, x^1, del, NULL) != 1) break;
get_available_cnt(g, x^1, del, &arc);
w = arc.v ^ 1;
if (get_available_cnt(g, w, del, NULL) != 1) break;
mark[x] = mark[w^1] = 1;
///l = asg_arc_len(arc_first(g, w));
get_available_cnt(g, w, del, &arc);
l = ((uint32_t)((arc).ul));
///w is the seq id + direction, l is the length of edge
///push element to the front of a queue
//kdq_unshift(uint64_t, q, (uint64_t)w<<32 | l);
q_occ++;
start = w, len += l;
x = w;
}
add_unitig:
if (start != UINT32_MAX) mark[start] = mark[end] = 1;
totalLen += len;
}
//kdq_destroy(uint64_t, q);
return totalLen;
}
void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen)
{
uint8_t *del = (uint8_t *)malloc(sizeof(uint8_t)*g->n_seq);
uint32_t v, k;
ma_utg_t* u = NULL;
memset(del, 1, g->n_seq);
(*primaryLen) = (*alterLen) = 0;
for (v = 0; v < ug->g->n_seq; ++v)
{
if(ug->g->seq[v].del) continue;
if(ug->g->seq[v].c == ALTER_LABLE) continue;
u = &(ug->u.a[v]);
if(u->m == 0) continue;
for (k = 0; k < u->n; k++)
{
del[u->a[k]>>33] = 0;
}
}
(*primaryLen) = get_specific_contig_length(g, del);
for (v = 0; v < g->n_seq; ++v)
{
del[v] = 1 - del[v];
}
(*alterLen) = get_specific_contig_length(g, del);
free(del);
}
int if_ploid_sample(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut,
hap_alignment_struct_pip* hap_buf, hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp,
uint32_t minLen, double purge_threshold)
{
asg_t* nsg = ug->g;
uint64_t v, k, total_bases = 0, alter_bases = 0, primary_bases = 0, purge_bases = 0;
kt_for(asm_opt.thread_num, hap_alignment_advance_worker, hap_buf, nsg->n_seq);
filter_hap_overlaps_by_length(all_ovlp, minLen);
normalize_hap_overlaps_advance(all_ovlp, back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
get_contig_length(ug, read_g, &primary_bases, &alter_bases);
total_bases = primary_bases + alter_bases;
fprintf(stderr, "primary_bases: %lu\n", primary_bases);
fprintf(stderr, "alter_bases: %lu\n", alter_bases);
fprintf(stderr, "total_bases: %lu\n", total_bases);
for (v = 0; v < all_ovlp->num; v++)
{
for (k = 0; k < all_ovlp->x[v].a.n; k++)
{
purge_bases += all_ovlp->x[v].a.a[k].x_end_pos - all_ovlp->x[v].a.a[k].x_beg_pos;
}
}
purge_bases = purge_bases/2;
fprintf(stderr, "purge_bases: %lu\n", purge_bases);
alter_bases = alter_bases + purge_bases;
fprintf(stderr, "new alter_bases: %lu\n", alter_bases);
for (v = 0; v < all_ovlp->num; v++)
{
all_ovlp->x[v].a.n = 0;
}
for (v = 0; v < back_all_ovlp->num; v++)
{
back_all_ovlp->x[v].a.n = 0;
}
if(alter_bases > total_bases * purge_threshold) return 1;
return 0;
}
void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density,
uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio,
@@ -3961,9 +4203,11 @@ uint32_t just_contain)
asg_arc_t* p = NULL;
int r;
hap_alignment_struct_pip hap_buf;
long long k_mer_only, coverage_only;
hap_buf.cov_threshold = get_read_coverage_thres(ug, read_g, ruIndex, position_index, sources, coverage_cut, read_g->n_seq, COV_COUNT);
fprintf(stderr, "cov_threshold: %lld\n", hap_buf.cov_threshold);
hap_buf.cov_threshold = get_read_coverage_thres(ug, read_g, ruIndex, position_index,
sources, coverage_cut, read_g->n_seq, COV_COUNT, &k_mer_only, &coverage_only);
///fprintf(stderr, "cov_threshold: %lld\n", hap_buf.cov_threshold);
@@ -3998,6 +4242,24 @@ uint32_t just_contain)
sources, reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp,
0.05, &all_ovlp);
if(hap_buf.cov_threshold < 0)
{
if(if_ploid_sample(ug, read_g, ruIndex, sources, reverse_sources, coverage_cut,
&hap_buf, &all_ovlp, &back_all_ovlp, purege_minLen, 0.333))
{
///if peak is het, coverage peak is more reliable
hap_buf.cov_threshold = coverage_only * HET_PEAK_RATE;
}
else
{
///if peak is homo, k-mer peak is more reliable
hap_buf.cov_threshold = k_mer_only * HOM_PEAK_RATE;
}
}
asm_opt.hom_global_coverage = hap_buf.cov_threshold;
fprintf(stderr, "cov_threshold: %lld\n", hap_buf.cov_threshold);
///kt_for(asm_opt.thread_num, hap_alignment_worker, &hap_buf, nsg->n_seq);
kt_for(asm_opt.thread_num, hap_alignment_advance_worker, &hap_buf, nsg->n_seq);
Purge_Dups.h
+5 -1
View File
@@ -7,6 +7,10 @@
#include "Overlaps.h"
#include "Hash_Table.h"
#define COV_COUNT 1024
#define HOM_PEAK_RATE 1.25
#define HET_PEAK_RATE (HOM_PEAK_RATE*2)
#define ALTER_COV_THRES 0.9
#define REAL_ALTER_THRES 0.1
void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density,
@@ -14,7 +18,7 @@ uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float
uint32_t just_contain);
void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge,
uint32_t is_circle, uint64_t* rLen);
void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen);
void enable_debug_mode(uint32_t mode);
#endif