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

better seeding

Browse Source
This commit is contained in:
chhylp123
2021-08-22 09:53:34 -04:00
parent bfff640a82
commit 37b07e4d33
15 changed files with 637 additions and 406 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+9 -6
View File
@@ -13,7 +13,7 @@
#include "rcut.h"
void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres,
int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct);
int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp);
void ha_get_ug_candidates(ha_abuf_t *ab, int64_t rid, ma_utg_t *u, ma_utg_v *ua, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag,
kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, double chain_match_rate);
void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list);
@@ -438,6 +438,7 @@ typedef struct {
kvec_t_u8_warp k_flag;
overlap_region tmp_region;
ma_utg_v *ua;
st_mt_t sp;
} ha_ovec_buf_t;
ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
@@ -454,6 +455,7 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
kv_init(b->b_buf.a);
kv_init(b->r_buf.a);
kv_init(b->k_flag.a);
kv_init(b->sp);
b->ab = ha_abuf_init();
if (!b->is_final) {
init_Cigar_record(&b->cigar1);
@@ -478,6 +480,7 @@ void ha_ovec_destroy(ha_ovec_buf_t *b)
kv_destroy(b->b_buf.a);
kv_destroy(b->r_buf.a);
kv_destroy(b->k_flag.a);
kv_destroy(b->sp);
if (!b->is_final) {
destory_Cigar_record(&b->cigar1);
destory_Graph(&b->POA_Graph);
@@ -541,7 +544,7 @@ static void worker_ovec(void *data, long i, int tid)
int fully_cov, abnormal;
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist,
0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
clear_Cigar_record(&b->cigar1);
clear_Round2_alignment(&b->round2);
@@ -599,7 +602,7 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
int fully_cov, abnormal, q_idx = k;
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist,
0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), &(R_INF_FLAG.candidate_count[q_idx]));
0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), &(R_INF_FLAG.candidate_count[q_idx]), &(b->sp));
clear_Cigar_record(&b->cigar1);
clear_Round2_alignment(&b->round2);
@@ -625,7 +628,7 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
pthread_mutex_lock(&R_INF_FLAG.OutputMutex);
fprintf(R_INF_FLAG.fp, "\n>%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i));
fprintf(R_INF_FLAG.fp, "\n>%.*s, rid: %ld\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i), i);
fprintf(R_INF_FLAG.fp, "%d-th round, len: %lu, hom_cov: %d, max_n_chain: %d\n",
asm_opt.number_of_round, Get_READ_LENGTH(R_INF, i), asm_opt.hom_cov, asm_opt.max_n_chain);
@@ -1291,7 +1294,7 @@ static void worker_ov_final(void *data, long i, int tid)
//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 0.001,
asm_opt.max_n_chain, 0, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
asm_opt.max_n_chain, 0, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
overlap_region_sort_y_id(b->olist.list, b->olist.length);
ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length);
@@ -1357,7 +1360,7 @@ static void worker_ov_final_high_het(void *data, long i, int tid)
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, HIGH_HET_ERROR_RATE,
asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
overlap_region_sort_y_id(b->olist.list, b->olist.length);
ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length);
CommandLines.cpp
+1
View File
@@ -154,6 +154,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->k_mer_length = 51;
asm_opt->hic_mer_length = 31;
asm_opt->mz_win = 51;
asm_opt->mz_rewin = 1000;
asm_opt->mz_sample_dist = 500;
asm_opt->bf_shift = 37;
asm_opt->max_kmer_cnt = 2000;
CommandLines.h
+2 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.15.5-r358"
#define HA_VERSION "0.15.5-r366"
#define VERBOSE 0
@@ -46,6 +46,7 @@ typedef struct {
int hic_mer_length;
int bub_mer_length;
int mz_win;
int mz_rewin;
int mz_sample_dist;
int bf_shift;
int max_kmer_cnt;
Hash_Table.cpp
+47 -5
View File
@@ -9,6 +9,8 @@ pthread_mutex_t output_mutex;
#define overlap_region_key(a) ((a).y_id)
KRADIX_SORT_INIT(overlap_region_sort, overlap_region, overlap_region_key, member_size(overlap_region, y_id))
#define normal_w(x, y) ((x)>=(y)?(x)/(y):1)
void overlap_region_sort_y_id(overlap_region *a, long long n)
{
radix_sort_overlap_region_sort(a, a + n);
@@ -467,7 +469,46 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp)
}
}
void print_chain(k_mer_hit* a, long long a_n, Chain_Data* dp, long long topN)
{
fprintf(stderr, "topN: %lld\n", topN);
long long max_score = -1, max_i = -1, max_n = 0;;
long long ss, i, j, current_j;
kvec_t(long long) si; kv_init(si);
for (ss = 0; ss < topN && ss < a_n; ss++){
for (i = 0, max_i = -1, max_score = -1; i < a_n; ++i) {
for (j = 0; j < (long long)si.n; j++){
if(i == si.a[j]) break;
}
if(j < (long long)si.n) continue;
if(dp->score[i] > max_score) max_score = dp->score[i], max_i = i;
}
if(max_i < 0) continue;
j = max_i; max_n = 0;
while (j >= 0)
{
current_j = j;
if(current_j == -1) continue;
j = dp->pre[j];
max_n++;
}
fprintf(stderr, "\nmax_i: %lld, max_score: %lld, max_n: %lld\n", max_i, max_score, max_n);
j = max_i;
while (j >= 0)
{
current_j = j;
if(current_j == -1) continue;
kv_push(long long, si, current_j);
j = dp->pre[j];
fprintf(stderr, "self_offset: %u, offset: %u, cnt: %u, score: %d\n",
a[current_j].self_offset, a[current_j].offset, a[current_j].cnt, dp->score[current_j]);
}
}
kv_destroy(si);
}
long long get_chainLen(long long x_beg, long long x_end, long long xLen,
@@ -556,7 +597,8 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc
break;
if (i < n_a) return -1;
bw_pen = 1.0 / bw_thres;
dp->score[0] = a[0].good? min_sc : min_sc>>1;
// dp->score[0] = a[0].good? min_sc : min_sc>>1;
dp->score[0] = normal_w(min_sc, (int64_t)a[0].cnt);
dp->pre[0] = -1, dp->indels[0] = 0, dp->self_length[0] = 0, dp->occ[0] = 1;
for (i = 1; i < n_a; ++i) {
int32_t score, dg;
@@ -570,7 +612,7 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc
dg = dx < dy? dx : dy;
if (dd > THRESHOLD_MAX_SIZE && dd > dg * bw_thres) break;
score = dg < min_sc? dg : min_sc;
if (!a[i].good) score >>= 1;
score = normal_w(score, (int64_t)a[i].cnt);
gap_rate = (double)tot_indel / tot_len;
score -= (int)(gap_rate * score * bw_pen);
dp->score[i] = dp->score[i-1] + score;
@@ -601,7 +643,6 @@ long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region*
resize_Chain_Data(dp, a_n);
ret = ha_chain_check(a, a_n, dp, min_score, band_width_threshold);
if (ret > 0) {
a_n = ret;
goto skip_dp;
@@ -617,7 +658,7 @@ long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region*
pos = a[i].offset;
self_pos = a[i].self_offset;
max_j = -1;
max_score = a[i].good? min_score : min_score>>1;
max_score = normal_w(min_score, (int64_t)a[i].cnt);
max_indels = 0;
max_self_length = 0;
@@ -648,7 +689,7 @@ long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region*
score = distance_min < min_score? distance_min : min_score;
///need to be fixed in r305
///if (!a[j].good) score = (score >> 1) + (score & 1);
if (!a[j].good) score >>= 1;
score = normal_w(score, (int64_t)a[j].cnt);
gap_rate = (double)((double)(total_indels)/(double)(total_self_length));
///if the gap rate > 0.06, score will be negative
@@ -684,6 +725,7 @@ long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region*
}
///debug_chain(a, a_n, dp);
// if((*result).x_id == 2162668 && (*result).y_id == 182804) print_chain(a, a_n, dp, 10);
skip_dp:
Hash_Table.h
+2 -2
View File
@@ -110,8 +110,8 @@ typedef struct
typedef struct
{
uint32_t readID:30, strand:1, good:1;
uint32_t offset, self_offset;
uint32_t readID:31, strand:1;
uint32_t offset, self_offset, cnt;
} k_mer_hit;
typedef struct {
Overlaps.cpp
+106 -224
View File
@@ -153,7 +153,7 @@ void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend)
///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)
R_to_U* ruIndex, uint8_t* is_r_het, 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;
@@ -187,7 +187,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
b_max.b_0 = b_0;
}
uint32_t max_count = 0, min_count = 0;
uint32_t max_count = 0, min_count = 0, n_het = 0, n_hom = 0;
ma_utg_t *node_min = NULL, *node_max = NULL;
if(ug != NULL)
{
@@ -217,7 +217,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
/************************BUG: don't forget****************************/
if(reverse_sources[qn].length > 0) min_count++;
///if(reverse_sources[qn].length >= 0) min_count++;
if((is_r_het[qn] & C_HET) || (is_r_het[qn] & P_HET)) n_het++;
n_hom++;
/************************BUG: don't forget****************************/
for (j = 0; j < (long long)reverse_sources[qn].length; j++)
{
@@ -268,7 +269,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
/************************BUG: don't forget****************************/
if(reverse_sources[qn].length > 0) min_count++;
///if(reverse_sources[qn].length >= 0) min_count++;
if((is_r_het[qn] & C_HET) || (is_r_het[qn] & P_HET)) n_het++;
n_hom++;
/************************BUG: don't forget****************************/
for (j = 0; j < (long long)reverse_sources[qn].length; j++)
@@ -300,7 +302,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
}
if(min_count == 0) return UNAVAILABLE;
if(max_count > min_count*asm_opt.purge_simi_thres/**DIFF_HAP_RATE**/) return PLOID;
if(max_count > min_count*asm_opt.purge_simi_thres && n_het >= n_hom*HET_HOM_RATE) return PLOID;
return NON_PLOID;
}
@@ -476,107 +478,6 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint32_t *min_count, uint32_t
(*max_count) = inp_match;
(*min_count) = hap_match;
}
inline uint32_t check_different_haps_base(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;
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, &tmp, &ELen_0, &max_stop_nodeLen, &max_stop_baseLen,
stops_threshold, b_0) == LOOP)
{
return UNAVAILABLE;
}
if(get_unitig(nsg, ug, v_1, &vEnd, &tmp, &ELen_1, &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_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****************************************/
calculate_match_cover(b_min.b_0->b.a, b_min.b_0->b.n, nsg, ug, read_sg,
reverse_sources, ruIndex, &min_count, &max_count);
/*****************************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****************************************/
calculate_match_cover(b_min.b_0->b.a, b_min.b_0->b.n, nsg, NULL, read_sg,
reverse_sources, ruIndex, &min_count, &max_count);
/*****************************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 == 67 && v_1 == 510)
// {
// fprintf(stderr, "v_0-%u, v_1-%u, min_count-%u, max_count-%u\n", v_0, v_1, min_count, max_count);
// }
if(min_count == 0) return UNAVAILABLE;
if(max_count > min_count*asm_opt.purge_simi_thres/**DIFF_HAP_RATE**/) return PLOID;
return NON_PLOID;
}
asg_t *asg_init(void)
{
@@ -12504,7 +12405,7 @@ ma_ug_t *ug, asg_t *read_sg, hap_cov_t *cov)
rid = (ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33));
uCov += cov->cov[rid];
if(t_ch) t_ch->is_r_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET;
if(t_ch) t_ch->ir_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET;
}
}
@@ -12520,7 +12421,7 @@ ma_ug_t *ug, asg_t *read_sg, hap_cov_t *cov)
rid = (ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33));
uLen += read_sg->seq[rid].len;
if(t_ch) t_ch->is_r_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET;
if(t_ch) t_ch->ir_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET;
}
if(occ >= thre_pri) break;
}
@@ -12693,7 +12594,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp)
void set_ug_coverage_aggressive(ma_ug_t *ug, uint32_t uID, asg_t* read_g,
const ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag,
trans_chain* t_ch, long long het_cov_thres)
uint8_t* is_r_het, long long het_cov_thres)
{
ma_utg_t *u = &(ug->u.a[uID]);
uint32_t k, j, rId, tn, is_Unitig;
@@ -12791,7 +12692,7 @@ trans_chain* t_ch, long long het_cov_thres)
if((R_bases <= 0) || ((C_bases/R_bases) <= het_cov_thres))
{
t_ch->is_r_het[rId] |= C_HET;
is_r_het[rId] |= C_HET;
}
}
@@ -12820,7 +12721,7 @@ trans_chain* t_ch, long long het_cov_thres)
}
void set_r_het_flag(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, trans_chain* t_ch)
void set_r_het_flag(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* is_r_het)
{
uint64_t m, dip_thre_max, dip_thres;
uint8_t* primary_flag = (uint8_t*)calloc(sg->n_seq, sizeof(uint8_t));
@@ -12842,7 +12743,7 @@ void set_r_het_flag(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, ma_hit_t_all
{
dip_thres = dip_thre_max;
///if(ug->u.a[m].n <= dip_thre_max) dip_thres = dip_thre_max * 1.1;
set_ug_coverage_aggressive(ug, m, sg, coverage_cut, sources, ruIndex, primary_flag, t_ch, dip_thres);
set_ug_coverage_aggressive(ug, m, sg, coverage_cut, sources, ruIndex, primary_flag, is_r_het, dip_thres);
}
free(primary_flag);
}
@@ -12859,7 +12760,6 @@ trans_chain* init_trans_chain(ma_ug_t *ug, uint64_t r_num)
memset(x->rUidx, -1, x->r_num*sizeof(uint32_t));
MALLOC(x->rUpos, r_num);
memset(x->rUpos, -1, x->r_num*sizeof(uint64_t));
CALLOC(x->is_r_het, x->r_num);
memset(&(x->b_buf_0), 0, sizeof(buf_t));
memset(&(x->b_buf_1), 0, sizeof(buf_t));
kv_init(x->topo_buf);
@@ -12939,7 +12839,6 @@ void destory_trans_chain(trans_chain **x)
kv_destroy((*x)->k_t_b);
free((*x)->rUidx);
free((*x)->rUpos);
free((*x)->is_r_het);
uint32_t k;
for (k = 0; k < (*x)->bed.n; k++) kv_destroy((*x)->bed.a[k]);
kv_destroy((*x)->bed);
@@ -13046,7 +12945,7 @@ void write_trans_chain(trans_chain* t_ch, const char *fn)
FILE* fp = fopen(buf, "w");
fwrite(&t_ch->r_num, sizeof(t_ch->r_num), 1, fp);
fwrite(t_ch->is_r_het, sizeof(uint8_t), t_ch->r_num, fp);
fwrite(t_ch->ir_het, sizeof(uint8_t), t_ch->r_num, fp);
uint32_t i;
fwrite(&t_ch->bed.n, sizeof(t_ch->bed.n), 1, fp);
@@ -13082,8 +12981,8 @@ trans_chain* load_hc_trans(const char *fn)
CALLOC(t_ch, 1);
flag += fread(&t_ch->r_num, sizeof(t_ch->r_num), 1, fp);
MALLOC(t_ch->is_r_het, t_ch->r_num);
flag += fread(t_ch->is_r_het, sizeof(uint8_t), t_ch->r_num, fp);
MALLOC(t_ch->ir_het, t_ch->r_num);
flag += fread(t_ch->ir_het, sizeof(uint8_t), t_ch->r_num, fp);
uint32_t i;
flag += fread(&t_ch->bed.n, sizeof(t_ch->bed.n), 1, fp);
@@ -13351,7 +13250,12 @@ int max_hang, int min_ovlp, R_to_U* ruIndex, bub_label_t* b_mask_t)
kv_init(new_rtg_edges.a); kv_init(d_edges.a);
ma_ug_t *ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
adjust_utg_advance(sg, ug, reverse_sources, ruIndex, b_mask_t);
uint8_t* is_r_het = NULL;
CALLOC(is_r_het, sg->n_seq);
set_r_het_flag(ug, sg, coverage_cut, sources, ruIndex, is_r_het);
adjust_utg_advance(sg, ug, reverse_sources, ruIndex, b_mask_t, is_r_het);
asg_t* nsg = (*ug).g;
uint32_t v, n_vtx = nsg->n_seq;
for (v = 0; v < n_vtx; ++v)
@@ -13366,6 +13270,7 @@ int max_hang, int min_ovlp, R_to_U* ruIndex, bub_label_t* b_mask_t)
kv_destroy(new_rtg_edges.a); kv_destroy(d_edges.a);
horder_clean_sg_by_utg(sg, ug);
free(is_r_het);
return ug;
}
@@ -13377,8 +13282,8 @@ int max_hang, int min_ovlp, R_to_U* ruIndex, bub_label_t* b_mask_t)
p->r_num = sg->n_seq; p->u_num = ug->u.n;
kv_malloc(p->bed, p->u_num); p->bed.n = p->u_num;
for (k = 0; k < p->bed.n; k++) kv_init(p->bed.a[k]);
CALLOC(p->is_r_het, p->r_num);
kv_u_trans_t *ta = get_utg_ovlp(&ug, sg, sources, reverse_sources, coverage_cut, ruIndex, max_hang, min_ovlp, NULL, b_mask_t, p->is_r_het);
CALLOC(p->ir_het, p->r_num);
kv_u_trans_t *ta = get_utg_ovlp(&ug, sg, sources, reverse_sources, coverage_cut, ruIndex, max_hang, min_ovlp, NULL, b_mask_t, p->ir_het);
p->k_trans = *ta; free(ta);
return p;
}
@@ -13539,7 +13444,7 @@ void set_trio_flag_by_cov(ma_ug_t *ug, asg_t *read_g, hap_cov_t *cov)
for (k = 0; k < u->n; k++)
{
if(R_INF.trio_flag[u->a[k]>>33]&SET_TRIO) continue;
if(cov->t_ch->is_r_het[u->a[k]>>33] == N_HET) continue;
if(cov->t_ch->ir_het[u->a[k]>>33] == N_HET) continue;
R_INF.trio_flag[u->a[k]>>33] |= flag;
}
}
@@ -13639,23 +13544,6 @@ const ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t*
return R_bases == 0? 0 : C_bases/R_bases;
}
void print_r_het(hap_cov_t *cov, uint8_t* trio_flag, const char* cmd)
{
if(cov && cov->t_ch)
{
fprintf(stderr, "\n+%s-is_r_het[1369536]=%u\n", cmd, cov->t_ch->is_r_het[1369536]);
fprintf(stderr, "+%s-is_r_het[5097804]=%u\n", cmd, cov->t_ch->is_r_het[5097804]);
fprintf(stderr, "+%s-is_r_het[603738]=%u\n", cmd, cov->t_ch->is_r_het[603738]);
}
if(trio_flag)
{
fprintf(stderr, "+%s-trio_flag[1369536]=%u\n", cmd, trio_flag[1369536]);
fprintf(stderr, "+%s-trio_flag[5097804]=%u\n", cmd, trio_flag[5097804]);
fprintf(stderr, "+%s-trio_flag[603738]=%u\n", cmd, trio_flag[603738]);
}
}
void kt_u_trans_t_idx(kv_u_trans_t *ta, uint32_t n)
{
radix_sort_u_trans(ta->a, ta->a + ta->n);
@@ -14774,7 +14662,7 @@ uint32_t cal_trio_vec(buf_t* b, ma_ug_t *ug, float thres)
}
int cut_trio_tip_primary(asg_t *g, ma_ug_t *ug, uint32_t max_ext, uint32_t trio_flag, uint32_t keep_out_node,
asg_t *read_sg, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint32_t min_edge_length)
asg_t *read_sg, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t* is_r_het, uint32_t min_edge_length)
{
double startTime = Get_T();
uint32_t n_vtx = g->n_seq * 2, v, w, i, cnt = 0, tipEvaluateLen, flag, inner_flag, operation, tip_trio_flag;
@@ -14856,8 +14744,8 @@ asg_t *read_sg, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint32_t min_e
if(operation == CUT) break;
if(aw[i].del) continue;
if(aw[i].v == (b.b.a[b.b.n-1]^1)) continue;
inner_flag = /**check_different_haps_base**/check_different_haps(g, ug, read_sg, b.b.a[b.b.n-1]^1, aw[i].v,
reverse_sources, &b_0, &b_1, ruIndex, min_edge_length, 1);
inner_flag = check_different_haps(g, ug, read_sg, b.b.a[b.b.n-1]^1, aw[i].v,
reverse_sources, &b_0, &b_1, ruIndex, is_r_het, min_edge_length, 1);
if(inner_flag == NON_PLOID) operation = CUT;
}
}
@@ -14989,8 +14877,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov, utg
n_reduced++;
operation = TRIM;
flag = /**check_different_haps_base**/check_different_haps(g, ug, read_sg, av[v_maxLen_i].v, av[i].v,
reverse_sources, &b_0, &b_1, ruIndex, min_edge_length, 1);
flag = check_different_haps(g, ug, read_sg, av[v_maxLen_i].v, av[i].v, reverse_sources,
&b_0, &b_1, ruIndex, cov->is_r_het, min_edge_length, 1);
// #define UNAVAILABLE (uint32_t)-1
// #define PLOID 0
// #define NON_PLOID 1
@@ -15096,8 +14984,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre
{
n_reduced++;
operation = TRIM;
flag = /**check_different_haps_base**/check_different_haps(g, ug, read_sg, a_convex[convex_i].v, a_convex[i].v,
reverse_sources, &b_0, &b_1, ruIndex, min_edge_length, stops_threshold);
flag = check_different_haps(g, ug, read_sg, a_convex[convex_i].v, a_convex[i].v,
reverse_sources, &b_0, &b_1, ruIndex, cov->is_r_het, min_edge_length, stops_threshold);
// #define UNAVAILABLE (uint32_t)-1
// #define PLOID 0
// #define NON_PLOID 1
@@ -15290,8 +15178,8 @@ hap_cov_t *cov, utg_trans_t *o)
if(return_flag != END_TIPS) continue;
flag = /**check_different_haps_base**/check_different_haps(g, ug, read_sg, av[base_maxLen_i].v, av[i].v,
reverse_sources, &b_0, &b_1, ruIndex, miniedgeLen, 1);
flag = check_different_haps(g, ug, read_sg, av[base_maxLen_i].v, av[i].v,
reverse_sources, &b_0, &b_1, ruIndex, cov->is_r_het, miniedgeLen, 1);
// if((av[i].v>>1) == 255 && (av[base_maxLen_i].v>>1) == 33)
// if((av[i].v>>1) == 1852 && (av[base_maxLen_i].v>>1) == 2441)
@@ -15616,10 +15504,8 @@ hap_cov_t *cov, utg_trans_t *o)
if(ll>convexLen && max_stop_baseLen>=ll*MAX_STOP_RATE)
{
// flag = check_different_haps(g, ug, read_sg, a_convex[convex_i].v, a_convex[i].v,
// reverse_sources, &b_0, &b_1, ruIndex, miniedgeLen, stops_threshold);
flag = /**check_different_haps_base**/check_different_haps(g, ug, read_sg, a_convex[convex_i].v, a_convex[i].v,
reverse_sources, &b_0, &b_1, ruIndex, miniedgeLen, stops_threshold);
flag = check_different_haps(g, ug, read_sg, a_convex[convex_i].v, a_convex[i].v,
reverse_sources, &b_0, &b_1, ruIndex, cov->is_r_het, miniedgeLen, stops_threshold);
// #define UNAVAILABLE (uint32_t)-1
// #define PLOID 0
// #define NON_PLOID 1
@@ -15682,7 +15568,7 @@ hap_cov_t *cov, utg_trans_t *o)
int detect_chimeric_by_topo(asg_t *g, ma_ug_t *ug, asg_t *read_sg,
ma_hit_t_alloc* reverse_sources, long long miniedgeLen, uint32_t stops_threshold, float drop_rate,
R_to_U* ruIndex, utg_trans_t *o)
R_to_U* ruIndex, utg_trans_t *o, uint8_t* is_r_het)
{
double startTime = Get_T();
uint32_t i, k, v_i, v_beg, v_end, selfLen, w1, w2, wv, nw, n_vtx = g->n_seq * 2, n_reduced = 0, convex, convex_T, read_num;
@@ -15762,8 +15648,8 @@ R_to_U* ruIndex, utg_trans_t *o)
}
if(k != b_0.b.n) break;
if(/**check_different_haps_base**/check_different_haps(g, ug, read_sg, wv, aw[i].v, reverse_sources, &b_0, &b_1,
ruIndex, miniedgeLen, stops_threshold)==PLOID)
if(check_different_haps(g, ug, read_sg, wv, aw[i].v, reverse_sources, &b_0, &b_1,
ruIndex, is_r_het, miniedgeLen, stops_threshold)==PLOID)
{
break;
}
@@ -15808,8 +15694,8 @@ R_to_U* ruIndex, utg_trans_t *o)
}
if(k != b_0.b.n) break;
if(/**check_different_haps_base**/check_different_haps(g, ug, read_sg, wv, aw[i].v, reverse_sources, &b_0, &b_1,
ruIndex, miniedgeLen, stops_threshold)==PLOID)
if(check_different_haps(g, ug, read_sg, wv, aw[i].v, reverse_sources, &b_0, &b_1,
ruIndex, is_r_het, miniedgeLen, stops_threshold)==PLOID)
{
break;
}
@@ -16055,7 +15941,7 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
/**********debug**********/
if(just_bubble_pop == 0)
{
cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, 2);
cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
}
/**********debug**********/
long long pre_cons = get_graph_statistic(g);
@@ -16073,7 +15959,7 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, cov, NULL);
asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL);
asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL);
detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL);
detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL, cov->is_r_het);
///need consider tangles
///note we need both the read graph and the untig graph
}
@@ -16082,14 +15968,14 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
}
if(just_bubble_pop == 0)
{
cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, 2);
cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
}
///print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, trio_flag, drop_ratio);
drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, DOUBLE_CHECK_THRES, asm_opt.trio_flag_occ_thres);
resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, trio_flag, drop_ratio);
drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het);
all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, cov->is_r_het, DOUBLE_CHECK_THRES, asm_opt.trio_flag_occ_thres);
if(is_first)
{
is_first = 0;
@@ -16133,24 +16019,23 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov)
asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, NULL, 1);
if(just_bubble_pop == 0)
{
cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2);
cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
}
// print_debug_gfa(read_g, ug, coverage_cut, "debug_init", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
long long pre_cons = get_graph_statistic(g);
long long cur_cons = 0;
while(pre_cons != cur_cons)
{
{
pre_cons = get_graph_statistic(g);
asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, NULL, 1);
if(just_bubble_pop == 0)
{
///need consider tangles
asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, NULL);
asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, NULL);
asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov, NULL);
asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL);
asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL);
detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL);
detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL, cov->is_r_het);
if(round != T_ROUND)
{
unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip,
@@ -16158,14 +16043,13 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov)
}
}
cur_cons = get_graph_statistic(g);
}
}
if(just_bubble_pop == 0)
{
cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex,
2);
cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
}
resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, (uint32_t)-1, drop_ratio);
drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, (uint32_t)-1, drop_ratio);
drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het);
// print_debug_gfa(read_g, ug, coverage_cut, "debug_clean_end", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, reverse_sources, 0, 1);
///print_graph_statistic(g, "end");
@@ -16196,7 +16080,7 @@ int min_ovlp, hap_cov_t *cov)
redo:
asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, o, 1);
cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2);
cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
long long pre_cons = get_graph_statistic(g);
long long cur_cons = 0;
@@ -16214,7 +16098,7 @@ int min_ovlp, hap_cov_t *cov)
asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov, o);
asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, o);
asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, o);
detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, o);
detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, o, cov->is_r_het);
cur_cons = get_graph_statistic(g);
}
@@ -16236,11 +16120,10 @@ int min_ovlp, hap_cov_t *cov)
}
cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex,
2);
cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2);
resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, (uint32_t)-1, drop_ratio);
drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, (uint32_t)-1, drop_ratio);
drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het);
print_debug_gfa(read_g, ug, coverage_cut, "debug_clean_end", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, reverse_sources, 0, 1);
@@ -16278,7 +16161,7 @@ void set_drop_trio_flag(ma_ug_t *ug)
}
void update_unitig_graph(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,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t* is_r_het,
uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate)
{
asg_t* nsg = ug->g;
@@ -16286,7 +16169,7 @@ uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate)
ma_utg_t *u;
uint8_t* primary_flag = (uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t));
drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex);
drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex, is_r_het);
while (n_reduce)
{
@@ -16365,7 +16248,7 @@ uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate)
}
}
drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex);
drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex, is_r_het);
asg_cleanup(nsg);
free(primary_flag);
}
@@ -16496,9 +16379,9 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex)
}
uint32_t unitig_simi(uint32_t x, uint32_t y, ma_ug_t* ug, ma_hit_t_alloc* reverse_sources,
R_to_U* ruIndex)
R_to_U* ruIndex, uint8_t* is_r_het)
{
uint32_t k, j, uId, tn, is_Unitig, rId, ref_unitig, min_count, max_count;
uint32_t k, j, uId, tn, is_Unitig, rId, ref_unitig, min_count, max_count, n_het, n_hom;
ma_utg_t *nsu_x = NULL, *nsu_y = NULL, *nsu_query = NULL;
nsu_x = &(ug->u.a[x]);
nsu_y = &(ug->u.a[y]);
@@ -16516,11 +16399,13 @@ R_to_U* ruIndex)
ref_unitig = x;
}
min_count = max_count = 0;
min_count = max_count = n_het = n_hom = 0;
for (k = 0; k < nsu_query->n; k++)
{
rId = nsu_query->a[k]>>33;
if(reverse_sources[rId].length >= 0) min_count++;
if((is_r_het[rId] & C_HET) || (is_r_het[rId] & P_HET)) n_het++;
n_hom++;
for (j = 0; j < reverse_sources[rId].length; j++)
{
@@ -16539,7 +16424,7 @@ R_to_U* ruIndex)
}
if(min_count == 0) return UNAVAILABLE;
if(max_count > min_count*asm_opt.purge_simi_thres/**DIFF_HAP_RATE**/) return PLOID;
if(max_count > min_count*asm_opt.purge_simi_thres && n_het >= n_hom*HET_HOM_RATE) return PLOID;
return NON_PLOID;
}
@@ -16573,7 +16458,7 @@ uint32_t* non_require, uint32_t* ambigious)
///note: to use this function, don't renew unitig graph!!!!!!!!!
void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate,
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate, int non_tig_occ)
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t* is_r_het, float double_check_rate, int non_tig_occ)
{
@@ -16679,7 +16564,7 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate, int n
if(flag_occ <= ((non_flag_occ+flag_occ)*drop_rate)) continue;
if((flag_occ+non_flag_occ) == 0) continue;
}**/
if(unitig_simi(uId, (uint32_t)(u_vecs.a.a[k]), ug, reverse_sources, ruIndex)==PLOID)
if(unitig_simi(uId, (uint32_t)(u_vecs.a.a[k]), ug, reverse_sources, ruIndex, is_r_het)==PLOID)
{
break;
}
@@ -16810,7 +16695,7 @@ ma_hit_t_alloc* 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 drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t* is_r_het)
{
uint32_t v, n_vtx = nsg->n_seq*2, convex_f, convex_b, i, nv;
long long ll, tmp, max_stop_nodeLen, max_stop_baseLen;
@@ -16845,8 +16730,8 @@ void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* re
}
get_real_length(nsg, convex_f, &convex_f);
if(convex_f != convex_b) continue;
if(/**check_different_haps_base**/check_different_haps(nsg, ug, read_g, v^1, av[i].v,
reverse_sources, &b_0, &b_1, ruIndex, 2, 1) == PLOID)
if(check_different_haps(nsg, ug, read_g, v^1, av[i].v,
reverse_sources, &b_0, &b_1, ruIndex, is_r_het, 2, 1) == PLOID)
{
av[i].del = 1;
asg_arc_del(nsg, av[i].v^1, v^1, 1);
@@ -16987,10 +16872,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
uint32_t v, n_vtx = nsg->n_seq;
hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources,
coverage_cut, max_hang, min_ovlp, asm_opt.purge_level_trio>0?1:0);
if(cov->t_ch)
{
set_r_het_flag(*ug, read_g, coverage_cut, sources, ruIndex, cov->t_ch);
}
if(cov->t_ch) cov->t_ch->ir_het = cov->is_r_het;
if(asm_opt.recover_atg_cov_min == -1024)
{
@@ -17011,9 +16893,9 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
}
adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t);
adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t, cov->is_r_het);
///primary_flag = get_utg_attributes(*ug, read_g, coverage_cut, sources, ruIndex);
update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0,
update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, cov->is_r_het, 0,
DOUBLE_CHECK_THRES, flag, drop_rate);
nsg = (*ug)->g;
@@ -17033,7 +16915,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
update_hap_label(*ug, read_g);
update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0,
update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, cov->is_r_het, 0,
DOUBLE_CHECK_THRES, flag, drop_rate);
force_trio_clean((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, flag, 0.55, 0.01, 5);
@@ -17059,7 +16941,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
///if(flag == MOTHER) print_untig_by_read(*ug, "m64043_200627_000137/124716590/ccs", 2789716, NULL, NULL, "beg");
update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 1,
update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, cov->is_r_het, 1,
FINAL_DOUBLE_CHECK_THRES, flag, drop_rate);
update_hap_label(NULL, read_g);
@@ -17730,7 +17612,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
ori = t_ch->b_buf_0.b.a[i]&1;
for (k = 0; k < p->n; k++)
{
t_ch->is_r_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
t_ch->ir_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
}
}
@@ -17744,7 +17626,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
ori = t_ch->b_buf_1.b.a[i]&1;
for (k = 0; k < p->n; k++)
{
t_ch->is_r_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
t_ch->ir_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
}
}
@@ -17776,7 +17658,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
ori = cov->t_ch->topo_res.a[i]&1;
for (k = 0; k < p->n; k++)
{
t_ch->is_r_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
t_ch->ir_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
}
/***********************x***********************/
@@ -17789,7 +17671,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
ori = t_ch->b_buf_0.b.a[k_i]&1;
for (k = 0; k < p->n; k++)
{
t_ch->is_r_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
t_ch->ir_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET;
}
}
/***********************y***********************/
@@ -20382,7 +20264,7 @@ uint32_t type)
inline uint32_t walk_through(asg_t *read_g, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, long long minLongUntig,
long long maxShortUntig, float l_untig_rate, float max_node_threshold, buf_t* b_0, buf_t* b_1,
kvec_t_u32_warp* u_vecs, uint8_t* visit, uint32_t v, uint32_t* r_beg, uint32_t* r_end,
uint32_t* r_next_uID, R_to_U* ruIndex)
uint32_t* r_next_uID, R_to_U* ruIndex, uint8_t* is_r_het)
{
(*r_beg) = (*r_end) = (uint32_t)-1;
asg_t* nsg = ug->g;
@@ -20478,8 +20360,8 @@ uint32_t* r_next_uID, R_to_U* ruIndex)
// #define UNAVAILABLE (uint32_t)-1
// #define PLOID 0
// #define NON_PLOID 1
if(returnFlag == 1 && /**check_different_haps_base**/check_different_haps(nsg, ug, read_g, beg, next_uID,
reverse_sources, b_0, b_1, ruIndex, minLongUntig-1, 1) == PLOID)
if(returnFlag == 1 && check_different_haps(nsg, ug, read_g, beg, next_uID,
reverse_sources, b_0, b_1, ruIndex, is_r_het, minLongUntig-1, 1) == PLOID)
{
///output_tangles(beg, next_uID, u_vecs->a.a, u_vecs->a.n, (char*)("???"));
returnFlag = 0;
@@ -22858,8 +22740,8 @@ float drop_ratio)
}
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)
long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint8_t* is_r_het,
uint32_t trio_flag, float drop_ratio)
{
buf_t b_0, b_1;
memset(&b_0, 0, sizeof(buf_t));
@@ -22910,7 +22792,7 @@ float drop_ratio)
{
flag = walk_through(read_g, ug, reverse_sources, minLongUntig,
maxShortUntig, l_untig_rate, max_node_threshold, &b_0, &b_1,
&u_vecs, visit, sv, &beg, &end, &next_uID, ruIndex);
&u_vecs, visit, sv, &beg, &end, &next_uID, ruIndex, is_r_het);
n_reduce += flag;
if(flag != UNROLL_M)
{
@@ -23604,7 +23486,7 @@ int get_arc_t(Edge_iter* x, asg_arc_t* get)
}
void unroll_simple_case_advance(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t, double dupLenThres)
void unroll_simple_case_advance(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t, uint8_t* is_r_het, double dupLenThres)
{
asg_t* nsg = ug->g;
uint32_t v, n_vtx = nsg->n_seq * 2, rnw, nw, beg, end, i;
@@ -23752,8 +23634,8 @@ void unroll_simple_case_advance(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* reve
continue;
}
if(/**check_different_haps_base**/check_different_haps(nsg, ug, read_g, beg, end, reverse_sources, &b_0, &b_1,
ruIndex, 2, 1) == PLOID)
if(check_different_haps(nsg, ug, read_g, beg, end, reverse_sources, &b_0, &b_1,
ruIndex, is_r_het, 2, 1) == PLOID)
{
continue;
}
@@ -23806,12 +23688,12 @@ void unroll_simple_case_advance(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* reve
void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t)
void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t, uint8_t* is_r_het)
{
double startTime = Get_T();
asg_t* nsg = ug->g;
unroll_simple_case_advance(ug, sg, reverse_sources, ruIndex, b_mask_t, 2.5);
drop_semi_circle(ug, ug->g, sg, reverse_sources, ruIndex);
unroll_simple_case_advance(ug, sg, reverse_sources, ruIndex, b_mask_t, is_r_het, 2.5);
drop_semi_circle(ug, ug->g, sg, reverse_sources, ruIndex, is_r_het);
asg_cleanup(nsg);
asg_symm(nsg);
///debug_utg_graph(ug, sg, 0, 0);
@@ -24455,6 +24337,8 @@ uint32_t is_collect_trans)
}
if(set) free(set);
CALLOC(x->is_r_het, read_g->n_seq);
set_r_het_flag(ug, read_g, coverage_cut, sources, ruIndex, x->is_r_het);
x->t_ch = NULL;
if(is_collect_trans) x->t_ch = init_trans_chain(ug, read_g->n_seq);
@@ -24468,6 +24352,7 @@ void destory_hap_cov_t(hap_cov_t **x)
{
free((*x)->cov);
free((*x)->pos_idx);
free((*x)->is_r_het);
kv_destroy((*x)->u_buffer.a);
kv_destroy((*x)->tailIndex.a);
kv_destroy((*x)->prevIndex.a);
@@ -24514,7 +24399,7 @@ void reset_trans_chain(trans_chain* t_ch, ma_utg_t *u)
{
uint32_t k = 0;
if(u->n == 0 || u->m == 0) return;
for (k = 0; k < u->n; k++) t_ch->is_r_het[u->a[k]>>33] = N_HET;
for (k = 0; k < u->n; k++) t_ch->ir_het[u->a[k]>>33] = N_HET;
}
void append_utg(ma_ug_t* ptg, ma_ug_t* atg, trans_chain* t_ch)
@@ -24733,8 +24618,8 @@ uint32_t collect_p_trans, uint32_t collect_p_trans_f)
ma_utg_t* u = NULL;
hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources,
coverage_cut, max_hang, min_ovlp, (asm_opt.purge_level_primary>0||i_cov)?1:0);
if(cov->t_ch) set_r_het_flag(*ug, read_g, coverage_cut, sources, ruIndex, cov->t_ch);
adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t);
if(cov->t_ch) cov->t_ch->ir_het = cov->is_r_het;
adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t, cov->is_r_het);
nsg = (*ug)->g;
n_vtx = nsg->n_seq;
@@ -28170,7 +28055,7 @@ void reset_bub(bubble_type* bub, ma_ug_t *ug, trans_chain* back_ug_chain, kvec_a
new_rtg_edges->a.n = 0;
///classify_untigs(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, max_hang, min_ovlp);
identify_bubbles(ug, bub, back_ug_chain->is_r_het, NULL);
identify_bubbles(ug, bub, back_ug_chain->ir_het, NULL);
// update_bubble_chain(ug, bub, 0, 1);
// resolve_bubble_chain_tangle(ug, bub);
// fprintf(stderr, "bub.f_bub: %lu, bub.b_bub: %lu, bub.b_end_bub: %lu, bub.tangle_bub: %lu, bub.cross_bub: %lu\n",
@@ -28461,7 +28346,7 @@ int max_hang, int min_ovlp, bubble_type* bub, long long gap_fuzz)
}
uint8_t *rescue_bubble_by_chain(asg_t *sg, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
void rescue_bubble_by_chain(asg_t *sg, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex,
float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, uint32_t chainLenThres, long long gap_fuzz,
bub_label_t* b_mask_t)
@@ -28505,14 +28390,11 @@ bub_label_t* b_mask_t)
rescue_missing_hap_ovlp(ug, sg, sources, coverage_cut, max_hang, min_ovlp, &bub, gap_fuzz);
}
uint8_t *het_flag = cov->t_ch->is_r_het;
cov->t_ch->is_r_het = NULL;
destory_bubbles(&bub);
destory_hap_cov_t(&cov);
ma_ug_destroy(ug);
kv_destroy(new_rtg_edges.a);
ma_ug_destroy(copy_ug); copy_ug = NULL;
return het_flag;
}
void update_unitig(long long step, long long init, ma_utg_t* nsu, asg_t *r_g,
@@ -30978,12 +30860,12 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex,
max_hang_length, mini_overlap_length);
ruIndex->is_het = rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3,
rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3,
ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz, &b_mask_t);
output_unitig_graph(sg, coverage_cut, o_file, sources, ruIndex, max_hang_length, mini_overlap_length);
// flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL;
flat_soma_v(sg, sources, ruIndex); free(ruIndex->is_het); ruIndex->is_het = NULL;
flat_soma_v(sg, sources, ruIndex);
output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources,
asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length);
Overlaps.h
+8 -7
View File
@@ -715,6 +715,7 @@ uint32_t stops_threshold, buf_t* b)
#define TRIO_DROP_LENGTH_THRES 0.8
#define MAX_STOP_RATE 0.6
#define TANGLE_MISSED_THRES 0.6
#define HET_HOM_RATE 0.7
typedef struct {
uint32_t father_occ;
@@ -747,17 +748,17 @@ typedef struct{
} bub_label_t;
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, bub_label_t* b_mask_t);
long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint8_t* is_r_het,
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, bub_label_t* b_mask_t, uint8_t* is_r_het);
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, 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, bub_label_t* b_mask_t);
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, uint32_t is_bubble_check, uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t);
void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate, int non_tig_occ);
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);
ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t* is_r_het, float double_check_rate, int non_tig_occ);
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, uint8_t* is_r_het);
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, bub_label_t* b_mask_t);
void get_unitig_trio_flag(ma_utg_t* nsu, uint32_t flag, uint32_t* require, uint32_t* non_require, uint32_t* ambigious);
@@ -822,7 +823,7 @@ typedef struct {
typedef struct{
uint32_t* rUidx;
uint64_t* rUpos;
uint8_t* is_r_het;
uint8_t* ir_het;
uint32_t r_num, u_num;
kvec_t(bed_in) bed;
kvec_t(uint32_t) topo_buf;
@@ -848,7 +849,7 @@ typedef struct {
kvec_asg_arc_t_offset u_buffer;
kvec_t_i32_warp tailIndex;
kvec_t_i32_warp prevIndex;
///hc_links* link;
uint8_t* is_r_het;
trans_chain* t_ch;
}hap_cov_t;
Purge_Dups.cpp
+13 -20
View File
@@ -176,6 +176,7 @@ long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* cove
///if the second peak is not significant
if(max2 < max * 0.05 || min > max2 * 0.95) max2 = max2_i = -1;
}
if(max2 < max*0.0075) max2 = max2_i = -1;
// look for smaller peak on the high end
@@ -205,7 +206,7 @@ long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* cove
if (max3 < max * 0.05 || min > max3 * 0.95 || max3_i > max_i * 3) max3 = max3_i = -1;
}
if (max3 < max*0.0075) max3 = max3_i = -1;
@@ -216,6 +217,13 @@ long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* cove
topo_peak = cov_buf[topo_peak_i];
if (topo_peak <= max * 0.05) topo_peak_i = topo_peak = -1;
}
if(asm_opt.purge_level_primary == 0)
{
(*hom_peak) = max_i;
return;
}
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;
@@ -345,12 +353,6 @@ ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
}
}
// if(qn == 1893151 || qn == 1929038)
// {
// fprintf(stderr, "qn: %lu, C_bases_primary: %lld, C_bases_alter: %lld, C_bases: %lld\n",
// qn, C_bases_primary, C_bases_alter, C_bases);
// }
C_bases = C_bases_primary + C_bases_alter;
if(C_bases_alter < C_bases * ALTER_COV_THRES) continue;
@@ -369,8 +371,6 @@ ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
}
}
///fprintf(stderr, "alter max_i: %lld, max: %lld\n", max_i, max);
///if(max_i < 5) max_i = max = -1;
if(alter_bases < primary_bases * REAL_ALTER_THRES) max_i = max = -1;
free(cov_buf);
@@ -4981,9 +4981,9 @@ void debug_p_g_t(p_g_t* pg, hap_cov_t *cov, asg_t *read_g)
offset += (uint32_t)u->a[i];
if(i >= sid && i <= eid)
{
if(cov->t_ch->is_r_het[u->a[i]>>33] != t->h_status)
if(cov->t_ch->ir_het[u->a[i]>>33] != t->h_status)
{
fprintf(stderr, "ERROR-(-3): is_r_het: %u, h_status: %u\n", cov->t_ch->is_r_het[u->a[i]>>33], t->h_status);
fprintf(stderr, "ERROR-(-3): is_r_het: %u, h_status: %u\n", cov->t_ch->ir_het[u->a[i]>>33], t->h_status);
}
}
}
@@ -5065,11 +5065,11 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k)
{
///if (k == u->n || (!!cov->t_ch->is_r_het[u->a[k]>>33]) != (!!cov->t_ch->is_r_het[u->a[l]>>33]))
if (k == u->n || cov->t_ch->is_r_het[u->a[k]>>33] != cov->t_ch->is_r_het[u->a[l]>>33])
if (k == u->n || cov->t_ch->ir_het[u->a[k]>>33] != cov->t_ch->ir_het[u->a[l]>>33])
{
kv_pushp(p_node_t, pg->pg_het_node, &t);
t->c_ug_id = uId;
t->h_status = cov->t_ch->is_r_het[u->a[l]>>33];
t->h_status = cov->t_ch->ir_het[u->a[l]>>33];
t->baseBeg = l_pos;
t->baseEnd = offset + read_g->seq[u->a[k-1]>>33].len - 1;
t->nodeBeg = l;
@@ -5261,7 +5261,6 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle
init_hap_overlaps_list(&all_ovlp, nsg->n_seq);
hap_overlaps_list back_all_ovlp;
init_hap_overlaps_list(&back_all_ovlp, nsg->n_seq);
///uint32_t junk_cov, hap_cov, dip_cov, junk_occ, repeat_occ, single_cov;
asg_arc_t t, *p = NULL;
int r;
hap_alignment_struct_pip hap_buf;
@@ -5382,12 +5381,6 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle
r = get_hap_arch(&(all_ovlp.x[uId].a.a[i]), ug->u.a[all_ovlp.x[uId].a.a[i].xUid].len,
ug->u.a[all_ovlp.x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
// if(all_ovlp.x[uId].a.a[i].xUid == 118 && all_ovlp.x[uId].a.a[i].yUid == 82)
// {
// fprintf(stderr, "r: %d\n", r);
// print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
// }
if(r < 0) continue;
p = asg_arc_pushp(pg->pg_h_lev);
*p = t;
Purge_Dups.h
+1 -1
View File
@@ -10,7 +10,7 @@
#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
#define REAL_ALTER_THRES 0.25
#define CHAIN_FILTER_RATE 0.7
#define SELF_EXIST 0
anchor.cpp
+25 -15
View File
@@ -1,4 +1,5 @@
#include <stdio.h>
#include <math.h>
#include "htab.h"
#include "ksort.h"
#include "Hash_Table.h"
@@ -7,8 +8,9 @@
typedef struct { // this struct is not strictly necessary; we can use k_mer_pos instead, with modifications
uint64_t srt;
uint32_t self_off:31, good:1;
uint32_t self_off;
uint32_t other_off;
uint32_t cnt;
} anchor1_t;
#define an_key1(a) ((a).srt)
@@ -23,7 +25,7 @@ KSORT_INIT(or_xs, overlap_region, oreg_xs_lt)
KSORT_INIT(or_ss, overlap_region, oreg_ss_lt)
typedef struct {
int n, good;
int n;
const ha_idxpos_t *a;
} seed1_t;
@@ -58,12 +60,13 @@ int ha_ov_type(const overlap_region *r, uint32_t len)
}
void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain,
kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct)
kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp)
{
uint32_t i, rlen;
uint64_t k, l;
double low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
double high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
if(low_occ < 2) low_occ = 2;
// prepare
clear_Candidates_list(cl);
@@ -73,7 +76,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
rlen = Get_READ_LENGTH(R_INF, rid); // read length
// get the list of anchors
ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, NULL);
ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin);
// minimizer of queried read
if (ab->mz.m > ab->old_mz_m) {
ab->old_mz_m = ab->mz.m;
@@ -83,7 +86,6 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
int n;
ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n);
ab->seed[i].n = n;
ab->seed[i].good = (n > low_occ && n < high_occ);
ab->n_a += n;
}
if (ab->n_a > ab->m_a) {
@@ -102,7 +104,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
uint8_t rev = z->rev == y->rev? 0 : 1;
an->other_off = y->pos;
an->self_off = rev? ucr->length - 1 - (z->pos + 1 - z->span) : z->pos;
an->good = s->good;
an->cnt = s->n;
an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off;
}
}
@@ -129,7 +131,16 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
p->strand = ab->a[k].srt >> 32 & 1;
p->offset = ab->a[k].other_off;
p->self_offset = ab->a[k].self_off;
p->good = ab->a[k].good;
if(ab->a[k].cnt > low_occ && ab->a[k].cnt < high_occ){
p->cnt = 1;
}
else if(ab->a[k].cnt <= low_occ){
p->cnt = 2;
}
else{
p->cnt = 1 + ((ab->a[k].cnt + (high_occ<<1) - 1)/(high_occ<<1));
p->cnt = pow(p->cnt, 1.1);
}
}
cl->length = ab->n_a;
@@ -269,7 +280,6 @@ kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f
int n;
ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n);
ab->seed[i].n = n;
ab->seed[i].good = 0;
ab->n_a += n;
}
if (ab->n_a > ab->m_a) {
@@ -288,7 +298,7 @@ kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f
uint8_t rev = z->rev == y->rev? 0 : 1;
an->other_off = y->pos;
an->self_off = rev? u->len - 1 - (z->pos + 1 - z->span) : z->pos;
an->good = s->good;
an->cnt = 1;
an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off;
}
}
@@ -315,7 +325,7 @@ kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f
p->strand = ab->a[k].srt >> 32 & 1;
p->offset = ab->a[k].other_off;
p->self_offset = ab->a[k].self_off;
p->good = ab->a[k].good;
p->cnt = 1;
}
cl->length = ab->n_a;
@@ -391,14 +401,14 @@ void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf)
void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres,
int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar,
kvec_t_u64_warp* dbg_ct)
kvec_t_u64_warp* dbg_ct, st_mt_t *sp)
{
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
extern void *ha_flt_tab_hp;
extern ha_pt_t *ha_idx_hp;
ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar, dbg_ct);
ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar, dbg_ct, sp);
if(ha_idx_hp)
{
@@ -428,7 +438,7 @@ void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overl
overlap_list->length = k;
ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab_hp, ha_idx_hp, f_cigar, dbg_ct);
ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab_hp, ha_idx_hp, f_cigar, dbg_ct, sp);
if(overlap_list->length + overlap_list_hp->length > overlap_list->size)
{
hic.cpp
+2 -2
View File
@@ -14786,7 +14786,7 @@ void debug_gfa_space(ha_ug_index* idx, ma_ug_t* ug, trans_chain* t_ch, kv_u_tran
memset(&bub, 0, sizeof(bubble_type));
bub.round_id = 0; bub.n_round = 2;
identify_bubbles(ug, &bub, t_ch->is_r_het, ref);
identify_bubbles(ug, &bub, t_ch->ir_het, ref);
hc_links link;
init_hc_links(&link, ug->g->n_seq, t_ch);
@@ -15687,7 +15687,7 @@ void resolve_tangles_hic(ha_ug_index *idx, bubble_type *bub, kvec_pe_hit *hits,
pe_hit *h_a = NULL;
u_trans_t *p = NULL;
identify_bubbles(idx->ug, bub, idx->t_ch->is_r_het, &(idx->t_ch->k_trans));
identify_bubbles(idx->ug, bub, idx->t_ch->ir_het, &(idx->t_ch->k_trans));
ta->idx.n = ta->n = 0;
if(hits->idx.n == 0) idx_hc_links(hits, idx, NULL);
htab.cpp
+4 -5
View File
@@ -675,7 +675,7 @@ static void worker_for_mz(void *data, long i, int tid)
ha_mz1_v *b = &s->mz_buf[tid];
s->mz_buf[tid].n = 0;
ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, asm_opt.mz_sample_dist, 0, 0,
(s->p->pt&&(s->p->flag&HAF_COUNT_REFINE))?s->p->pt:NULL, s->p->opt->min_rcnt, asm_opt.dp_min_len, asm_opt.dp_e, s->mt?&(s->mt[tid]):NULL);
(s->p->pt&&(s->p->flag&HAF_COUNT_REFINE))?s->p->pt:NULL, s->p->opt->min_rcnt, asm_opt.dp_min_len, asm_opt.dp_e, &(s->mt[tid]), asm_opt.mz_rewin);
s->mz[i].n = s->mz[i].m = b->n;
MALLOC(s->mz[i].a, b->n);
memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t));
@@ -807,17 +807,16 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
// s->mz && s->mz_buf are lists of minimzer vectors
CALLOC(s->mz, s->n_seq);
CALLOC(s->mz_buf, p->opt->n_thread);
s->mt = NULL;
if(s->p->pt&&(s->p->flag&HAF_COUNT_REFINE)) CALLOC(s->mt, p->opt->n_thread);
CALLOC(s->mt, p->opt->n_thread);
///calculate minimzers for each read, each read corresponds to one thread
kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq);
for (i = 0; i < p->opt->n_thread; ++i)
{
if(s->mt) free(s->mt[i].a);
free(s->mt[i].a);
free(s->mz_buf[i].a);
}
if(s->mt) free(s->mt);
free(s->mt);
free(s->mz_buf);
// insert minimizers
if (p->pt && !(p->flag&HAF_COUNT_REFINE)) {///insert whole minimizer
htab.h
+1 -1
View File
@@ -101,7 +101,7 @@ double yak_cpu_usage(void);
void ha_triobin(const hifiasm_opt_t *opt);
void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt);
void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws);
int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het);
void print_hist_lines(int n_cnt, int start_cnt, const int64_t *cnt);
void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs);
rcut.cpp
+6 -6
View File
@@ -889,9 +889,9 @@ trans_chain* t_ch)
offset += (uint32_t)u->a[i];
if(i >= sid && i <= eid)
{
if(t_ch->is_r_het[u->a[i]>>33] != t->hs)
if(t_ch->ir_het[u->a[i]>>33] != t->hs)
{
fprintf(stderr, "ERROR-d: is_r_het: %u, h_status: %u\n", t_ch->is_r_het[u->a[i]>>33], t->hs);
fprintf(stderr, "ERROR-d: is_r_het: %u, h_status: %u\n", t_ch->ir_het[u->a[i]>>33], t->hs);
}
}
}
@@ -939,11 +939,11 @@ void update_mc_edges(mc_g_t *mg, hap_overlaps_list* ha, kv_u_trans_t *ta, trans_
u = &(mg->ug->u.a[v]);
for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k)
{
if (k == u->n || t_ch->is_r_het[u->a[k]>>33] != t_ch->is_r_het[u->a[l]>>33])
if (k == u->n || t_ch->ir_het[u->a[k]>>33] != t_ch->ir_het[u->a[l]>>33])
{
kv_pushp(mc_interval_t, p, &t);
t->uID = v;
t->hs = t_ch->is_r_het[u->a[l]>>33];
t->hs = t_ch->ir_het[u->a[l]>>33];
t->bS = l_pos;
t->bE = offset + mg->rg->seq[u->a[k-1]>>33].len - 1;
@@ -2680,10 +2680,10 @@ void set_p_flag(mc_g_t *mg, uint32_t uID, uint8_t* trio_flag, trans_chain* t_ch,
for (i = 0; i < u->n; i++)
{
trio_flag[u->a[i]>>33] |= SET_TRIO;
if(t_ch->is_r_het[u->a[i]>>33] == N_HET) continue;
if(t_ch->ir_het[u->a[i]>>33] == N_HET) continue;
if(s == 0)
{
if(t_ch->is_r_het[u->a[i]>>33]&P_HET)//special case
if(t_ch->ir_het[u->a[i]>>33]&P_HET)//special case
{
trio_flag[u->a[i]>>33] |= FATHER;
}
sketch.cpp
+410 -111
View File
@@ -229,7 +229,7 @@ inline void hf_select(ha_mz1_v *p, int32_t si, int32_t ei, int32_t n, int32_t le
p->a[b[j].pos].rid = 0;
}
static void select_mz(ha_mz1_v *p, int len, int sample_dist, int32_t dp_min_len)
void select_mz(ha_mz1_v *p, int len, int sample_dist, int32_t dp_min_len)
{ // for high-occ minimizers, choose up to max_high_occ in each high-occ streak
int32_t i, last0 = -1, n = (int32_t)p->n, m = 0, nw[2], min_len;
ha_mz1_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation
@@ -289,6 +289,288 @@ static void select_mz(ha_mz1_v *p, int len, int sample_dist, int32_t dp_min_len)
p->n = n;
}
static inline int mzcmp_l(const ha_mz1_v *p, int32_t ai, int32_t bi)
{
if(ai >= 0 && bi >= 0){
ha_mz1_t *a = &(p->a[ai]), *b = &(p->a[bi]);
if(a->rid > 0 && b->rid > 0) return mzcmp(a, b);
return (a->rid == 0) - (b->rid == 0);
}
return (ai < 0) - (bi < 0);
}
#define GL(x, i) ((int64_t)((uint32_t)((x).a[(i)])))
#define A_M(p, i) ((i) >= 0 && (p).a[(i)].rid > 0)
int32_t qfw(ha_mz1_v *p, st_mt_t *mt, int32_t n, int32_t tot_l, int32_t ws, int32_t i, int32_t *mi)
{
int32_t m, si;
for (si = i, (*mi) = -1; i < n; i++){
if(GL(*mt, i) >= ws || (i+1 < n && GL(*mt, i) < ws && GL(*mt, i+1) > ws) ||
(i+1 == n && tot_l >= ws && GL(*mt, i) < ws)){
for (m = si; m <= i; m++){
if(!A_M(*p, m)) continue;
if(mzcmp_l(p, *mi, m) >= 0) (*mi) = m;
}
if((*mi) >= 0 && A_M(*p, *mi)){
for (m = si; m <= i; m++){
if(!A_M(*p, m)) continue;
if(mzcmp_l(p, *mi, m) == 0) mt->a[m] |= 0x100000000;
}
}
break;
}
}
return i;
}
void dbg_boundary(ha_mz1_v *p, st_mt_t *mt, int32_t w, int32_t k, int32_t tot_l)
{
if(tot_l < w + k -1) return;
int32_t i, m, n = p->n, s, a;
for (i = 0; i < n; i++){
if(GL(*mt, i) >= w+k-1){
for (m = s = a = 0; m <= i; m++){
if(!A_M(*p, m)) continue;
if(GL(*mt, m) <= w+k-1){
a++;
if(mt->a[m]&0x100000000) s++;
}
}
if(a > 0 && s == 0){
fprintf(stderr, "\nERROR1, s: %d, n: %d, tot_l: %d, end_l: %ld\n", s, n, tot_l, GL(*mt, i));
for (m = s = a = 0; m <= i; m++){
if(!A_M(*p, m)) continue;
if(GL(*mt, m) <= w+k-1){
fprintf(stderr, "lp: %ld\n", GL(*mt, m));
a++;
if(mt->a[m]&0x100000000) s++;
}
}
}
break;
}
}
if(i == n){
for (m = s = a = 0; m < n; m++){
if(!A_M(*p, m)) continue;
if(GL(*mt, m) <= w+k-1){
a++;
if(mt->a[m]&0x100000000) s++;
}
}
if(a > 0 && s == 0) fprintf(stderr, "ERROR2\n");
}
for (i = n-1; i >= 0; i--)
{
if (GL(*mt, i) + w <= tot_l + 1) {
for (m = i, s = a = 0; m < n; m++){
if(!A_M(*p, m)) continue;
if(GL(*mt, m) + w >= tot_l + 1){
a++;
if(mt->a[m]&0x100000000) s++;
}
}
if(a > 0 && s == 0) {
fprintf(stderr, "\nERROR3, s: %d, n: %d, tot_l: %d, end_l: %ld\n", s, n, tot_l, GL(*mt, i));
for (m = i, s = a = 0; m < n; m++){
if(!A_M(*p, m)) continue;
if(GL(*mt, m) + w >= tot_l + 1){
fprintf(stderr, "lp: %ld\n", GL(*mt, m));
a++;
if(mt->a[m]&0x100000000) s++;
}
}
}
break;
}
}
if(i < 0){
for (m = s = a = 0; m < n; m++){
if(!A_M(*p, m)) continue;
if(GL(*mt, m) + w >= tot_l + 1){
a++;
if(mt->a[m]&0x100000000) s++;
}
}
if(a > 0 && s == 0) fprintf(stderr, "ERROR4\n");
}
}
static void select_mz_h(ha_mz1_v *p, st_mt_t *mt, int len, int sample_dist, int32_t w, int32_t k, int32_t tot_l)
{ // for high-occ minimizers, choose up to max_high_occ in each high-occ streak
int32_t i, mi = -1, si, last0 = -1, n = (int32_t)p->n, m = 0, ws = w + k - 1;
if (n == 0) return;
assert(n < 1<<27); // 27 is the number of bits for ha_mz1_t::pos; this should be safe as there are more bases than minimizers
for (i = m = 0, last0 = -1; i <= n; ++i) {
if (i == n || p->a[i].rid == 0) {
if (i - last0 > 1) {
int32_t ps = last0 < 0? 0 : p->a[last0].pos;
int32_t pe = i == n? len : p->a[i].pos;
if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){
last0 = -2;
m++;
break;
}
}
last0 = i;
}
}
if (m == 0) return; // no high-frequency k-mers; do nothing
if(last0 >= -1) goto ff;
i = 0;
i = qfw(p, mt, n, tot_l, ws, i, &mi);
if(i == n) goto ff;
for (si = 0, i++; i < n; i++){
for (; si < i; si++){
if(GL(*mt, si) + w > GL(*mt, i)) break;
}
// a new minimum; then write the old min
if(mzcmp_l(p, i, mi) <= 0) {
if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
mi = i;
}// old min has moved outside the window
else if(si > mi){
if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
for (m = si, mi = -1; m <= i; m++){
if(mzcmp_l(p, mi, m) >= 0) mi = m;
}
if(A_M(*p, mi)){
for (m = si; m <= i; m++){
if(!A_M(*p, m)) continue;
if(mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000;
}
}
}
}
if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
for (i = n - 1; si < n && GL(*mt, si) + w <= tot_l + 1; si++){
if(si > mi){
if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
for (m = si, mi = -1; m <= i; m++){
if(mzcmp_l(p, mi, m) >= 0) mi = m;
}
if(A_M(*p, mi)){
for (m = si; m <= i; m++){
if(!A_M(*p, m)) continue;
if(mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000;
}
}
}
}
/**
dbg_boundary(p, mt, w, k, tot_l);
fprintf(stderr, "\n");
for (i = 0; i < (int32_t)p->n; ++i){
if(p->a[i].rid == 0) continue;
fprintf(stderr, "%cl: %u, pos: %lu, cnt: %lu, key: %lu, i: %d\n", "+-"[!!(mt->a[i]&0x100000000)],
(uint32_t)mt->a[i], p->a[i].pos, p->a[i].rid, p->a[i].x, i);
// if (mt->a[i]&0x100000000){
// fprintf(stderr, "+l: %u, pos: %lu, cnt: %lu\n", (uint32_t)mt->a[i], p->a[i].pos, p->a[i].rid);
// }
}
**/
ha_mz1_t b[MAX_MAX_HIGH_OCC];
for (i = 0, last0 = -1; i <= n; ++i) {
if (i == n || p->a[i].rid == 0) {
if (i - last0 > 1) {
int32_t ps = last0 < 0? 0 : p->a[last0].pos;
int32_t pe = i == n? len : p->a[i].pos;
if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){
for (m = last0 + 1, mi = 0; m < i; ++m){
if(mt->a[m]&0x100000000) p->a[m].rid = 0, mi++;
}
if(mi == 0) hf_select(p, last0, i, n, len, sample_dist, b, 0);
}
}
last0 = i;
}
}
ff:
for (i = n = 0; i < (int32_t)p->n; ++i) // squeeze out filtered minimizers
if (p->a[i].rid == 0)
p->a[n++] = p->a[i];
p->n = n;
}
void debug_pl(const char *str, int len, int w, int k, int is_hpc, ha_mz1_v *p, const void *hf, st_mt_t *mt)
{
int i, l, dbi, dbcnt = 0, kmer_span = 0;
tiny_queue_t tq;
memset(&tq, 0, sizeof(tiny_queue_t));
uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};
for (i = l = dbi = 0; i < len; ++i) {
int c = seq_nt4_table[(uint8_t)str[i]];
if (c < 4) { // not an ambiguous base
int z;
if (is_hpc) {
int skip_len = 1;
if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {
for (skip_len = 2; i + skip_len < len; ++skip_len)
if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)
break;
i += skip_len - 1; // put $i at the end of the current homopolymer run
}
tq_push(&tq, skip_len);
kmer_span += skip_len;
///how many bases that are covered by this HPC k-mer
///kmer_span includes at most k HPC elements
if (tq.count > k) kmer_span -= tq_shift(&tq);
} else kmer_span = l + 1 < k? l + 1 : k;
///kmer_span should be used for HPC k-mer
///non-HPC k-mer, kmer_span should be k
///kmer_span is used to calculate anchor pos on reverse complementary strand
kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer
kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer
kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1;
if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand
z = kmer[1] < kmer[3]? 0 : 1; // strand
++l;
if (l >= k && kmer_span < 256) {
uint64_t y;
int32_t cnt;
y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
cnt = hf? ha_ft_cnt(hf, y) : 0;
for (dbi = 0; dbi < mt->n; dbi++)
{
if(p->a[dbi].x == y && p->a[dbi].rid == cnt && p->a[dbi].pos == i && p->a[dbi].rev == z && p->a[dbi].span == kmer_span)
{
if(l != (int)mt->a[dbi]) fprintf(stderr, "ERROR\n");
dbcnt++;
}
}
}
} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
}
if(dbcnt != mt->n) fprintf(stderr, "ERROR\n");
if(mt->n != (int)p->n) fprintf(stderr, "ERROR\n");
for (dbi = 1; dbi < mt->n; dbi++)
{
if(p->a[dbi].pos <= p->a[dbi-1].pos || (int)mt->a[dbi] <= (int)mt->a[dbi-1])
{
fprintf(stderr, "ERROR\n");
}
}
}
/**
* Find symmetric (w,k)-minimizers on a DNA sequence
*
@@ -300,123 +582,140 @@ static void select_mz(ha_mz1_v *p, int len, int sample_dist, int32_t dp_min_len)
* @param is_hpc homopolymer-compressed or not
* @param p minimizers
*/
void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt)
{ ///in default, w = 51, k = 51, is_hpc = 1
/**
uint64_t x;
uint64_t rid:28, pos:27, rev:1, span:8;
**/
extern void *ha_ct_table;
static const ha_mz1_t dummy = { UINT64_MAX, (1<<28) - 1, 0, 0 };
uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};
int i, j, l, buf_pos, min_pos, kmer_span = 0;
ha_mz1_t buf[256], min = dummy;
tiny_queue_t tq;
void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws)
{ ///in default, w = 51, k = 51, is_hpc = 1
/**
uint64_t x;
uint64_t rid:28, pos:27, rev:1, span:8;
**/
extern void *ha_ct_table;
static const ha_mz1_t dummy = { UINT64_MAX, (1<<28) - 1, 0, 0, 0};
uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};
int i, j, l, tl = 0, buf_pos, min_pos, kmer_span = 0;
ha_mz1_t buf[256], min = dummy;
uint32_t buf_p[256], min_s = (uint32_t)-1;
tiny_queue_t tq;
assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
if (dbg_ct != NULL) dbg_ct->a.n = 0;
if (k_flag != NULL) {
kv_resize(uint8_t, k_flag->a, (uint64_t)len);
k_flag->a.n = len;
memset(k_flag->a.a, 0, k_flag->a.n);
}
assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
if (dbg_ct != NULL) dbg_ct->a.n = 0;
if (k_flag != NULL) {
kv_resize(uint8_t, k_flag->a, (uint64_t)len);
k_flag->a.n = len;
memset(k_flag->a.a, 0, k_flag->a.n);
}
memset(buf, 0xff, w * sizeof(ha_mz1_t));
memset(&tq, 0, sizeof(tiny_queue_t));
///len/w is the evaluated minimizer numbers
kv_resize(ha_mz1_t, *p, p->n + len/w);
memset(buf, 0xff, w * sizeof(ha_mz1_t));
memset(&tq, 0, sizeof(tiny_queue_t));
///len/w is the evaluated minimizer numbers
kv_resize(ha_mz1_t, *p, p->n + len/w);
kv_resize(uint64_t, *mt, (int64_t)p->m); mt->n = p->n;
for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
int c = seq_nt4_table[(uint8_t)str[i]];
ha_mz1_t info = dummy;
if (c < 4) { // not an ambiguous base
int z;
if (is_hpc) {
int skip_len = 1;
if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {
for (skip_len = 2; i + skip_len < len; ++skip_len)
if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)
break;
i += skip_len - 1; // put $i at the end of the current homopolymer run
}
tq_push(&tq, skip_len);
kmer_span += skip_len;
///how many bases that are covered by this HPC k-mer
///kmer_span includes at most k HPC elements
if (tq.count > k) kmer_span -= tq_shift(&tq);
} else kmer_span = l + 1 < k? l + 1 : k;
///kmer_span should be used for HPC k-mer
///non-HPC k-mer, kmer_span should be k
///kmer_span is used to calculate anchor pos on reverse complementary strand
for (i = l = tl = buf_pos = min_pos = 0; i < len; ++i) {
int c = seq_nt4_table[(uint8_t)str[i]];
ha_mz1_t info = dummy;
if (c < 4) { // not an ambiguous base
int z;
if (is_hpc) {
int skip_len = 1;
if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {
for (skip_len = 2; i + skip_len < len; ++skip_len)
if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)
break;
i += skip_len - 1; // put $i at the end of the current homopolymer run
}
tq_push(&tq, skip_len);
kmer_span += skip_len;
///how many bases that are covered by this HPC k-mer
///kmer_span includes at most k HPC elements
if (tq.count > k) kmer_span -= tq_shift(&tq);
} else kmer_span = l + 1 < k? l + 1 : k;
///kmer_span should be used for HPC k-mer
///non-HPC k-mer, kmer_span should be k
///kmer_span is used to calculate anchor pos on reverse complementary strand
if (k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC
if (k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC
kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer
kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer
kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1;
if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand
z = kmer[1] < kmer[3]? 0 : 1; // strand
++l;
if (l >= k && kmer_span < 256) {
uint64_t y;
int32_t cnt, filtered;
y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
cnt = hf? ha_ft_cnt(hf, y) : 0;
filtered = (cnt >= 1<<28);
if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));
if (!filtered) info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span; // initially ha_mz1_t::rid keeps the k-mer count
if (k_flag != NULL) k_flag->a.a[i]++;
if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++;
}
} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer
kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer
kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1;
if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand
z = kmer[1] < kmer[3]? 0 : 1; // strand
++l; tl++;
if (l >= k && kmer_span < 256) {
uint64_t y;
int32_t cnt, filtered;
y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
cnt = hf? ha_ft_cnt(hf, y) : 0;
filtered = (cnt >= 1<<28);
if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));
if (!filtered) info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span; // initially ha_mz1_t::rid keeps the k-mer count
if (k_flag != NULL) k_flag->a.a[i]++;
if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++;
}
} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
//for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i
//i is the real base iterator, while l is the HPC base iterator
//only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX)
//but even if l < k, infor is still stored into buf
buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
for (j = buf_pos + 1; j < w; ++j)
if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
for (j = 0; j < buf_pos; ++j)
if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
}
/**
* There are three cases:
* 1. info.x <= min.x, means info is a new minimizer
* 2. info.x > min.x, info is not a new minimizer
* (1) buf_pos != min_pos, do nothing
* (2) buf_pos == min_pos, means current minimizer has moved outside the window
* **/
///three cases: 1.
if (info.x <= min.x) { // a new minimum; then write the old min
if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
min = info, min_pos = buf_pos;
} else if (buf_pos == min_pos) { // old min has moved outside the window
if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
///buf_pos == min_pos, means current minimizer has moved outside the window
///so for now we need to find a new minimizer at the current window (w k-mers)
for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers
if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
for (j = 0; j <= buf_pos; ++j)
if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j;
buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
buf_p[buf_pos] = l;
if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
for (j = buf_pos + 1; j < w; ++j){
if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){
kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);
}
}
for (j = 0; j < buf_pos; ++j){
if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){
kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);
}
}
}
if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
for (j = 0; j <= buf_pos; ++j)
if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
}
}
if (++buf_pos == w) buf_pos = 0;
}
if (min.x != UINT64_MAX)
kv_push(ha_mz1_t, *p, min);
if (sample_dist > w) select_mz(p, len, MAX_HIGH_OCC, dp_min_len);
if (dp_min_len > 0 && pt && mt) refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt);
for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts
p->a[i].rid = rid;
/**
* There are three cases:
* 1. info.x <= min.x, means info is a new minimizer
* 2. info.x > min.x, info is not a new minimizer
* (1) buf_pos != min_pos, do nothing
* (2) buf_pos == min_pos, means current minimizer has moved outside the window
* **/
///three cases: 1.
if (mzcmp(&min, &info) >= 0) { // a new minimum; then write the old min
if (l >= w + k && min.x != UINT64_MAX){
kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s);
}
min = info, min_pos = buf_pos, min_s = buf_p[buf_pos];
} else if (buf_pos == min_pos) { // old min has moved outside the window
if (l >= w + k - 1 && min.x != UINT64_MAX){
kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s);
}
///buf_pos == min_pos, means current minimizer has moved outside the window
///so for now we need to find a new minimizer at the current window (w k-mers)
for (j = buf_pos + 1, min = dummy; j < w; ++j) // the two loops are necessary when there are identical k-mers
if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j, min_s = buf_p[j]; // >= is important s.t. min is always the closest k-mer
for (j = 0; j <= buf_pos; ++j)
if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j, min_s = buf_p[j];
if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){
kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);
}
for (j = 0; j <= buf_pos; ++j)
if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){
kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);
}
}
}
if (++buf_pos == w) buf_pos = 0;
}
if (min.x != UINT64_MAX){
kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s);
}
// debug_pl(str, len, w, k, is_hpc, p, hf, mt);
// if (sample_dist > w) select_mz(p, len, MAX_HIGH_OCC, dp_min_len);
select_mz_h(p, mt, len, sample_dist, ws, k, tl);
if (dp_min_len > 0 && pt && mt) refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt);
for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts
p->a[i].rid = rid;
}
void ha_sketch_worse(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt)