Compare commits

...
92 Commits
Author SHA1 Message Date
chhylp123 284cd0784a Merge branch 'master' into hifiasm-v0.14 2021-02-09 18:46:29 -05:00
chhylp123 ceeb4562af fix hc_edge 2021-02-09 18:44:06 -05:00
chhylp123 5bb22f701c debug hic_hits 2021-02-09 18:40:49 -05:00
chhylp123 8467d58ef1 for v0.14 2021-02-09 13:00:23 -05:00
chhylp123 fda13cb0e4 update -z 2021-02-02 12:45:50 -05:00
chhylp123 02b9a3e109 hic multiple rounds 2021-01-31 18:49:08 -05:00
chhylp123 75ce7cd576 for v0.14-hic 2021-01-31 10:34:21 -05:00
chhylp123 a009026e65 fix bug for broken bubbles 2021-01-15 07:33:13 -05:00
chhylp123 88f5b11f04 bubble graph 2021-01-01 04:50:42 -05:00
chhylp123 c545fd84c7 hic bubble 2020-12-29 09:15:44 -05:00
chhylp123 5d506e82c9 phasing 2020-12-28 03:35:17 -05:00
chhylp123 53a655d802 init partition 2020-12-26 19:01:24 -05:00
chhylp123 f3e390eee8 for hic 2020-12-17 02:53:31 -05:00
chhylp123 73b5ef6769 bubble phasing 2020-11-21 20:50:31 -05:00
chhylp123 3a1a3750e9 Merge branch 'master' of https://github.com/chhylp123/Long_read_assembly 2020-11-16 01:49:07 -05:00
chhylp123 d87314a07e resolve conflict 2020-11-16 01:44:12 -05:00
chhylp123 ab073aa1fc hic first 2020-11-14 23:35:32 -05:00
Heng Li eddb8173da added --min-hist-cnt for #49 2020-10-26 16:28:30 -04:00
chhylp123 35e3423f93 debug counting 2020-10-26 12:32:55 -04:00
chhylp123 0cc02035c5 HiC full index 2020-10-26 11:01:48 -04:00
chhylp123 b2b4624cea rerelease v0.13 2020-10-22 14:24:59 -04:00
chhylp123 de4163e881 Merge branch 'master' of https://github.com/chhylp123/Long_read_assembly
merge README
2020-10-21 22:03:43 -04:00
chhylp123 53adde2b37 fix conflict 2020-10-21 22:01:07 -04:00
chhylp123 91d4915da7 release v0.13 2020-10-21 21:57:43 -04:00
chhylp123 04f94514ee fix circle 2020-10-19 19:51:17 -04:00
chhylp123 b1f844801a fix c->l bug 2020-10-19 16:39:53 -04:00
chhylp123 e5d503ea65 update polishing 2020-10-19 13:22:45 -04:00
chhylp123 9d4e3b3283 update polishing 2020-10-18 17:59:31 -04:00
chhylp123 753ea81194 reomve redundancy 2020-10-16 18:08:12 -04:00
chhylp123 89fcb35a64 for debugging 2020-10-11 14:32:33 -04:00
Heng Li 8cb131dca3 added a "Why hifiasm?" section 2020-10-05 01:01:58 -04:00
Heng Li 00128c3076 Fixed the wrong download link 2020-10-04 23:01:16 -04:00
Heng Li 519d5b1eb1 Updated the download link 2020-10-04 23:00:27 -04:00
chhylp123 7a314d1ef3 fix type in README 2020-08-27 23:57:48 -04:00
chhylp123 55f95c0191 fix typo 2020-08-26 00:25:42 -04:00
chhylp123 b629dd8556 freeze-0.12-r304 2020-08-23 23:38:04 -04:00
chhylp123 0a2559e4a0 freeze-r304 2020-08-23 23:36:45 -04:00
chhylp123 724990cbe7 freeze-r304 2020-08-23 23:34:36 -04:00
chhylp123 03cf0d475c save bfore multiple overlaps 2020-08-23 08:46:33 -04:00
chhylp123 90d3cb69b9 fix version number 2020-08-05 18:04:47 -04:00
chhylp123 927d66262b rerelase v0.11 2020-08-05 17:50:53 -04:00
chhylp123 e031c70d0f Merge branch 'master' of https://github.com/chhylp123/Long_read_assembly 2020-08-05 17:46:17 -04:00
chhylp123 a8aa37d201 fix makefile 2020-08-05 17:45:23 -04:00
Heng Li ebedae2a21 updated Makefile 2020-08-05 17:33:14 -04:00
chhylp123 320dc58127 fix bug at r_utg 2020-08-04 16:43:50 -04:00
chhylp123 ef5ba0cc2c update help/man 2020-07-31 21:32:02 -04:00
chhylp123 8bdc21799c update READMEN 2020-07-31 12:42:27 -04:00
chhylp123 7c385e219a update version number 2020-07-31 12:36:37 -04:00
chhylp123 d92f9519ae Merge branch 'master' of https://github.com/chhylp123/Long_read_assembly 2020-07-31 12:34:49 -04:00
chhylp123 6ca9d01f06 update version number 2020-07-31 12:28:30 -04:00
chhylp123 23f0a8b020 Merge pull request #23 from chhylp123/hifiasm_high_het
Hifiasm high het
2020-07-31 12:24:43 -04:00
chhylp123 9d659fd2e0 update v0.10 2020-07-31 12:23:25 -04:00
chhylp123 ace34ae2b8 for high-het and coverage 2020-07-28 13:07:35 -04:00
Heng Li 94008c02ac Re-release hifiasm-0.9 (r289) 2020-07-19 14:44:28 -04:00
Heng Li 1e4d67ee64 r288: reorganized command line help messages 2020-07-19 13:43:30 -04:00
Heng Li ff4ef89fe3 r287: fixed the INT32_MAX compiling error 2020-07-19 11:22:23 -04:00
chhylp123 952471657f change max coverage of pri-range 2020-07-18 17:40:58 -04:00
chhylp123 88ab801717 update HA_VERSION 2020-07-13 22:45:49 -04:00
chhylp123 a309e610ca resolve conflict 2020-07-13 22:44:40 -04:00
chhylp123 fda858f9e7 update coverage 2020-07-13 22:34:26 -04:00
chhylp123 c01b593d7c using coverage 2020-07-12 16:03:01 -04:00
chhylp123 a25c37ed89 update k-mer peaks 2020-07-07 13:20:25 -04:00
chhylp123 2e56ef7db5 apply coverage to purge_dups 2020-07-06 10:02:03 -04:00
Heng Li 309e885365 Release hifiasm-0.8 (r279) 2020-06-27 19:56:49 -04:00
Heng Li 5cb42cafc2 more README changees 2020-06-27 19:25:38 -04:00
chhylp123 428cc01d10 update link 2020-06-27 15:33:49 -04:00
chhylp123 3481a2fe2b reorganize readme 2020-06-27 15:28:36 -04:00
Heng Li f1e78720aa more getting started; added NA12878 and links 2020-06-27 14:22:39 -04:00
chhylp123 e401f88e57 update test datasets 2020-06-27 13:03:46 -04:00
chhylp123 7a2e02cb77 fix typo 2020-06-27 12:56:02 -04:00
chhylp123 3a88079a6e update typo 2020-06-27 12:44:12 -04:00
chhylp123 be26acc0e0 update readme 2020-06-27 12:43:14 -04:00
chhylp123 03c9028848 update readme 2020-06-27 12:36:33 -04:00
chhylp123 0071251c82 update README 2020-06-27 12:35:47 -04:00
chhylp123 39cf3f5dd8 update README 2020-06-27 12:34:41 -04:00
chhylp123 acba806ce1 update readme 2020-06-27 12:30:57 -04:00
chhylp123 46987e792a update readme 2020-06-27 12:29:23 -04:00
chhylp123 1ec650e3e5 update readme 2020-06-27 12:27:16 -04:00
chhylp123 61657f4b29 update readme 2020-06-27 12:06:59 -04:00
chhylp123 35ee94f040 update README 2020-06-27 11:46:48 -04:00
chhylp123 ed137c2ac6 update readme 2020-06-27 11:38:08 -04:00
chhylp123 9cc97472d9 update readme 2020-06-27 11:25:21 -04:00
chhylp123 449950cddb update readme 2020-06-27 10:57:22 -04:00
chhylp123 71f6c93d58 add options 2020-06-08 17:51:36 -04:00
Heng Li 97f6efc4da r255: faster way to extract reads 2020-06-03 21:38:59 -04:00
Heng Li eddbc61952 r254: added --ex-list and --ex-iter 2020-06-03 20:46:34 -04:00
Heng Li e247c5f042 r253: Merge branch 'master' into dev-lh3 2020-05-23 21:43:22 -04:00
chhylp123 d36e1e7782 fixed bugs for trio-binning 2020-05-17 11:53:12 -04:00
Heng Li f6ea45b58a r251: always keep the best hits 2020-05-08 00:40:53 -04:00
Heng Li e544360043 r250: fixed a bug in selecting top chains 2020-05-07 14:49:49 -04:00
Heng Li 3364a29b87 r249: make max overlap error rate CLI options 2020-05-06 01:16:50 -04:00
Heng Li 331e437c72 r248: contained reads not handled properly 2020-05-03 00:27:53 -04:00
29 changed files with 21412 additions and 5432 deletions
+533 -74
View File
@@ -11,9 +11,12 @@
#include "htab.h"
#include "kthread.h"
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);
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);
void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list);
All_reads R_INF;
Debug_reads R_INF_FLAG;
void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int* new_length)
{
@@ -296,30 +299,6 @@ void push_overlaps(ma_hit_t_alloc* paf, overlap_region_alloc* overlap_list, int
}
}
int if_exact_match(char* x, long long xLen, char* y, long long yLen, long long xBeg, long long xEnd, long long yBeg, long long yEnd)
{
long long overlapLen = xEnd - xBeg + 1;
if(yEnd - yBeg + 1 == overlapLen)
{
long long i;
for (i = 0; i < overlapLen; i++)
{
if(x[xBeg + i] != y[yBeg + i])
{
break;
}
}
if(i == overlapLen)
{
return 1;
}
}
return 0;
}
long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, overlap_region_alloc* overlap_list, int flag)
{
@@ -343,6 +322,64 @@ long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_l
///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand
tmp.rev = overlap_list->list[i].y_pos_strand;
/**********************target***************************/
tmp.tn = overlap_list->list[i].y_id;
if(tmp.rev == 1)
{
long long y_readLen = R_INF.read_length[overlap_list->list[i].y_id];
tmp.ts = y_readLen - overlap_list->list[i].y_pos_e - 1;
tmp.te = y_readLen - overlap_list->list[i].y_pos_s - 1;
}
else
{
tmp.ts = overlap_list->list[i].y_pos_s;
tmp.te = overlap_list->list[i].y_pos_e;
}
///the end pos is open
tmp.te++;
/**********************target***************************/
tmp.bl = R_INF.read_length[overlap_list->list[i].y_id];
tmp.ml = overlap_list->list[i].strong;
tmp.no_l_indel = overlap_list->list[i].without_large_indel;
tmp.el = overlap_list->list[i].shared_seed;
add_ma_hit_t_alloc(paf, &tmp);
}
}
return available_overlaps;
}
long long push_final_overlaps_increment(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, overlap_region_alloc* overlap_list, int flag)
{
long long i = 0;
long long available_overlaps = paf->length;
ma_hit_t tmp;
///clear_ma_hit_t_alloc(paf); // paf has been preallocated, so we don't need preallocation
for (i = 0; i < (long long)overlap_list->length; i++)
{
if (overlap_list->list[i].is_match == flag)
{
available_overlaps++;
/**********************query***************************/
//the interval of overlap is half-open [start, end)
tmp.qns = overlap_list->list[i].x_id;
tmp.qns = tmp.qns << 32;
tmp.qns = tmp.qns | (uint64_t)(overlap_list->list[i].x_pos_s);
///the end pos is open
tmp.qe = overlap_list->list[i].x_pos_e + 1;
/**********************query***************************/
///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand
tmp.rev = overlap_list->list[i].y_pos_strand;
@@ -383,6 +420,7 @@ typedef struct {
UC_Read self_read, ovlp_read;
Candidates_list clist;
overlap_region_alloc olist;
overlap_region_alloc olist_hp;
ha_abuf_t *ab;
// error correction related buffers
int64_t num_read_base, num_correct_base, num_recorrect_base;
@@ -392,6 +430,10 @@ typedef struct {
Correct_dumy correct;
haplotype_evdience_alloc hap;
Round2_alignment round2;
kvec_t_u32_warp b_buf;
kvec_t_u64_warp r_buf;
kvec_t_u8_warp k_flag;
overlap_region tmp_region;
} ha_ovec_buf_t;
ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
@@ -403,6 +445,11 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
init_UC_Read(&b->ovlp_read);
init_Candidates_list(&b->clist);
init_overlap_region_alloc(&b->olist);
init_overlap_region_alloc(&b->olist_hp);
init_fake_cigar(&(b->tmp_region.f_cigar));
kv_init(b->b_buf.a);
kv_init(b->r_buf.a);
kv_init(b->k_flag.a);
b->ab = ha_abuf_init();
if (!b->is_final) {
init_Cigar_record(&b->cigar1);
@@ -421,7 +468,12 @@ void ha_ovec_destroy(ha_ovec_buf_t *b)
destory_UC_Read(&b->ovlp_read);
destory_Candidates_list(&b->clist);
destory_overlap_region_alloc(&b->olist);
destory_overlap_region_alloc(&b->olist_hp);
ha_abuf_destroy(b->ab);
destory_fake_cigar(&(b->tmp_region.f_cigar));
kv_destroy(b->b_buf.a);
kv_destroy(b->r_buf.a);
kv_destroy(b->k_flag.a);
if (!b->is_final) {
destory_Cigar_record(&b->cigar1);
destory_Graph(&b->POA_Graph);
@@ -451,6 +503,7 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b)
{
int64_t i, mem = 0, mem_clist, mem_olist;
mem_clist = b->clist.size * sizeof(k_mer_hit) + b->clist.chainDP.size * 7 * 4;
mem_olist = b->olist.size * sizeof(overlap_region);
for (i = 0; i < (int64_t)b->olist.size; ++i) {
const overlap_region *r = &b->olist.list[i];
@@ -458,6 +511,14 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b)
mem_olist += r->f_cigar.size * 8;
mem_olist += r->boundary_cigars.size * sizeof(window_list);
}
mem_olist += b->olist_hp.size * sizeof(overlap_region);
for (i = 0; i < (int64_t)b->olist_hp.size; ++i) {
const overlap_region *r = &b->olist_hp.list[i];
mem_olist += r->w_list_size * sizeof(window_list);
mem_olist += r->f_cigar.size * 8;
mem_olist += r->boundary_cigars.size * sizeof(window_list);
}
mem = ha_abuf_mem(b->ab) + mem_clist + mem_olist;
if (!b->is_final) {
mem += sizeof(Cigar_record) + b->cigar1.lost_base_size + b->cigar1.size * 4;
@@ -475,7 +536,8 @@ static void worker_ovec(void *data, long i, int tid)
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
int fully_cov, abnormal;
ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
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);
clear_Cigar_record(&b->cigar1);
clear_Round2_alignment(&b->round2);
@@ -497,36 +559,133 @@ static void worker_ovec(void *data, long i, int tid)
}
R_INF.paf[i].is_abnormal = abnormal;
if (b->save_ov) {
R_INF.trio_flag[i] = AMBIGU;
///need to be fixed in r305
// if(ha_idx_hp == NULL)
// {
// R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &(b->k_flag), RESEED_HP_RATE, Get_READ_LENGTH(R_INF, i), NULL);
// }
if (R_INF.trio_flag[i] != AMBIGU || b->save_ov) {
int is_rev = (asm_opt.number_of_round % 2 == 0);
push_overlaps(&(R_INF.paf[i]), &b->olist, 1, &R_INF, is_rev);
push_overlaps(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, is_rev);
}
}
static void worker_ovec_related_reads(void *data, long i, int tid)
{
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
int required_read_name_length = strlen(asm_opt.required_read_name);
uint64_t k;
if (required_read_name_length == (int)Get_NAME_LENGTH((R_INF),i)
&&
memcmp(asm_opt.required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
{
ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i));
recover_UC_Read(&b->self_read, &R_INF, i);
fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq);
uint64_t k, queryNameLen;
for (k = 0; k < R_INF_FLAG.query_num; k++)
{
queryNameLen = strlen(R_INF_FLAG.read_name[k]);
if (queryNameLen != Get_NAME_LENGTH((R_INF),i)) continue;
if (memcmp(R_INF_FLAG.read_name[k], Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
{
break;
}
}
for (k = 0; k < b->olist.length; k++) {
fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
recover_UC_Read(&b->self_read, &R_INF, b->olist.list[k].y_id);
fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq);
if(k < R_INF_FLAG.query_num)
{
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]));
clear_Cigar_record(&b->cigar1);
clear_Round2_alignment(&b->round2);
correct_overlap(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
&b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal);
b->num_read_base += b->self_read.length;
b->num_correct_base += b->correct.corrected_base;
b->num_recorrect_base += b->round2.dumy.corrected_base;
push_cigar(R_INF.cigars, i, &b->cigar1);
push_cigar(R_INF.second_round_cigar, i, &b->round2.cigar);
R_INF.paf[i].is_fully_corrected = 0;
if (fully_cov) {
if (get_cigar_errors(&b->cigar1) == 0 && get_cigar_errors(&b->round2.cigar) == 0)
R_INF.paf[i].is_fully_corrected = 1;
}
R_INF.paf[i].is_abnormal = abnormal;
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, "%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);
fprintf(R_INF_FLAG.fp, "***************************k-mer counts (%d)***************************\n", (int)(R_INF_FLAG.candidate_count[q_idx].a.n));
sort_kvec_t_u64_warp(&(R_INF_FLAG.candidate_count[q_idx]), 0);
for (k = 0; k < R_INF_FLAG.candidate_count[q_idx].a.n; k++)
{
fprintf(R_INF_FLAG.fp, "[%lu] Count(%u): %lu, filtered: %lu\n", k,
(uint32_t)R_INF_FLAG.candidate_count[q_idx].a.a[k], R_INF_FLAG.candidate_count[q_idx].a.a[k]>>33,
(R_INF_FLAG.candidate_count[q_idx].a.a[k]>>32)&(uint64_t)1);
}
fprintf(R_INF_FLAG.fp, "***************************forward ovlp***************************\n");
for (k = 0; k < b->olist.length; k++)
{
if(b->olist.list[k].is_match != 1) continue;
fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e,
b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
}
}
fprintf(R_INF_FLAG.fp, "***************************reverse ovlp***************************\n");
for (k = 0; k < b->olist.length; k++)
{
if(b->olist.list[k].is_match != 2) continue;
fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e,
b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
}
fprintf(R_INF_FLAG.fp, "***************************unmatched ovlp***************************\n");
for (k = 0; k < b->olist.length; k++)
{
if(b->olist.list[k].is_match == 1) continue;
if(b->olist.list[k].is_match == 2) continue;
fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e,
b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
}
R_INF.trio_flag[i] = AMBIGU;
///need to be fixed in r305
// if(ha_idx_hp == NULL)
// {
// R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &(b->k_flag), RESEED_HP_RATE, Get_READ_LENGTH(R_INF, i), R_INF_FLAG.fp);
// }
fprintf(R_INF_FLAG.fp, "R_INF.trio_flag[%ld]: %u\n", i, R_INF.trio_flag[i]);
pthread_mutex_unlock(&R_INF_FLAG.OutputMutex);
}
}
static inline long long get_N_occ(char* seq, long long length)
{
long long j, N_occ = 0;
@@ -630,25 +789,76 @@ void Output_corrected_reads()
fclose(output_file);
}
void debug_print_pob_regions()
{
uint64_t i, total = 0;
for (i = 0; i < R_INF.total_reads; i++)
{
if(R_INF.trio_flag[i]!=AMBIGU)
{
total++;
fprintf(stderr, "(%lu) %.*s\n", i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
}
}
fprintf(stderr, "total hp reads: %lu, R_INF.total_reads: %lu\n", total, R_INF.total_reads);
exit(1);
}
void rescue_hp_reads(ha_ovec_buf_t **b)
{
int hom_cov, het_cov;
ha_flt_tab_hp = ha_idx_hp = NULL;
if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
ha_flt_tab_hp = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 1);
}
ha_idx_hp = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 1, &R_INF, &hom_cov, &het_cov);
if (asm_opt.required_read_name)
kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads);
else
kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads);
ha_ft_destroy(ha_flt_tab_hp); ha_flt_tab_hp = NULL;
ha_pt_destroy(ha_idx_hp); ha_idx_hp = NULL;
}
void ha_overlap_and_correct(int round)
{
int i, hom_cov;
int i, hom_cov, het_cov, r_out = 0;
ha_ovec_buf_t **b;
ha_ecsave_buf_t *e;
ha_flt_tab_hp = ha_idx_hp = NULL;
if((ha_idx == NULL)&&(asm_opt.flag & HA_F_VERBOSE_GFA)&&(round == asm_opt.number_of_round - 1))
{
r_out = 1;
}
if(asm_opt.required_read_name) init_Debug_reads(&R_INF_FLAG, asm_opt.required_read_name); // for debugging only
// overlap and correct reads
CALLOC(b, asm_opt.thread_num);
for (i = 0; i < asm_opt.thread_num; ++i)
b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1));
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF, &hom_cov); // build the index
if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
if(ha_idx) hom_cov = asm_opt.hom_cov;
if(ha_idx == NULL) ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
///debug_adapter(&asm_opt, &R_INF);
if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
ha_opt_update_cov(&asm_opt, hom_cov);
if (asm_opt.required_read_name)
kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads);
else
kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads);
if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
ha_pt_destroy(ha_idx);
ha_idx = 0;
ha_idx = NULL;
// collect statistics
for (i = 0; i < asm_opt.thread_num; ++i) {
@@ -660,8 +870,8 @@ void ha_overlap_and_correct(int round)
}
free(b);
if (asm_opt.required_read_name) exit(0); // for debugging only
if (asm_opt.required_read_name) destory_Debug_reads(&R_INF_FLAG), exit(0); // for debugging only
// save corrected reads to R_INF
CALLOC(e, asm_opt.thread_num);
for (i = 0; i < asm_opt.thread_num; ++i) {
@@ -677,8 +887,10 @@ void ha_overlap_and_correct(int round)
free(e[i].second_round_read);
}
free(e);
///debug_print_pob_regions();
}
void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf,
UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact)
{
@@ -760,6 +972,116 @@ UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact)
}
}
int check_chain_indels(Fake_Cigar* chain, long long xBeg, long long xEnd, float indel_rate)
{
uint64_t i = 0;
long long indels = 0, xOffset;
if(chain->length != 0)
{
indels += abs(get_fake_gap_shift(chain, 0));
xOffset = get_fake_gap_pos(chain, 0);
if(indels > (xOffset - xBeg + 1) * indel_rate) return 0;
for (i = 1; i < chain->length; i++)
{
indels += abs((get_fake_gap_shift(chain, i) - get_fake_gap_shift(chain, i-1)));
xOffset = get_fake_gap_pos(chain, i);
if(indels > (xOffset - xBeg + 1) * indel_rate) return 0;
}
}
if(indels > (xEnd - xBeg + 1) * indel_rate) return 0;
return 1;
}
void update_overlaps_chain_width(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf,
UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact, float indel_rate)
{
uint64_t inner_j = 0;
uint64_t j = 0;
long long x_overlapLen, y_overlapLen;
while (j < overlap_list->length && inner_j < paf->length)
{
if(overlap_list->list[j].y_id < paf->buffer[inner_j].tn)
{
j++;
}
else if(overlap_list->list[j].y_id > paf->buffer[inner_j].tn)
{
inner_j++;
}
else
{
if(check_chain_indels(&(overlap_list->list[j].f_cigar), overlap_list->list[j].x_pos_s,
overlap_list->list[j].x_pos_e, indel_rate) == 1)
{
if(overlap_list->list[j].y_pos_strand == paf->buffer[inner_j].rev)
{
x_overlapLen = Get_qe(paf->buffer[inner_j]) - Get_qs(paf->buffer[inner_j]) + 1;
y_overlapLen = Get_te(paf->buffer[inner_j]) - Get_ts(paf->buffer[inner_j]) + 1;
if(x_overlapLen < y_overlapLen) x_overlapLen = y_overlapLen;
x_overlapLen = x_overlapLen * 0.1;
// if(
// ((DIFF(overlap_list->list[j].x_pos_s, Get_qs(paf->buffer[inner_j])) < x_overlapLen)
// && (DIFF(overlap_list->list[j].x_pos_e, Get_qe(paf->buffer[inner_j])) < x_overlapLen))
// ||
// ((DIFF(overlap_list->list[j].y_pos_s, Get_ts(paf->buffer[inner_j])) < x_overlapLen)
// && (DIFF(overlap_list->list[j].y_pos_e, Get_te(paf->buffer[inner_j])) < x_overlapLen)))
if(
((DIFF(overlap_list->list[j].x_pos_s, Get_qs(paf->buffer[inner_j])) < (uint64_t)x_overlapLen)
&& (DIFF(overlap_list->list[j].x_pos_e, Get_qe(paf->buffer[inner_j])) < (uint64_t)x_overlapLen))
||
((DIFF(overlap_list->list[j].y_pos_s, Get_ts(paf->buffer[inner_j])) < (uint64_t)x_overlapLen)
&& (DIFF(overlap_list->list[j].y_pos_e, Get_te(paf->buffer[inner_j])) < (uint64_t)x_overlapLen))
)
{
overlap_list->list[j].is_match = is_match;
overlap_list->list[j].strong = paf->buffer[inner_j].ml;
overlap_list->list[j].without_large_indel = paf->buffer[inner_j].no_l_indel;
if(is_exact == 1)
{
if(overlap_list->list[j].y_pos_strand == 0)
{
recover_UC_Read(overlap_read, &R_INF, overlap_list->list[j].y_id);
}
else
{
recover_UC_Read_RC(overlap_read, &R_INF, overlap_list->list[j].y_id);
}
if(if_exact_match(g_read->seq, g_read->length, overlap_read->seq, overlap_read->length,
overlap_list->list[j].x_pos_s, overlap_list->list[j].x_pos_e,
overlap_list->list[j].y_pos_s, overlap_list->list[j].y_pos_e))
{
overlap_list->list[j].shared_seed = 1;
}
else
{
overlap_list->list[j].shared_seed = 0;
}
}
}
else
{
overlap_list->list[j].is_match = 3;
}
}
else
{
overlap_list->list[j].is_match = 3;
}
}
j++;
inner_j++;
}
}
}
void update_exact_overlaps(overlap_region_alloc* overlap_list, UC_Read* g_read, UC_Read* overlap_read)
{
uint64_t j;
@@ -767,6 +1089,12 @@ void update_exact_overlaps(overlap_region_alloc* overlap_list, UC_Read* g_read,
{
if (overlap_list->list[j].is_match != 1)
{
if((overlap_list->list[j].x_pos_e + 1 - overlap_list->list[j].x_pos_s) !=
(overlap_list->list[j].y_pos_e + 1 - overlap_list->list[j].y_pos_s))
{
continue;
}
if(overlap_list->list[j].y_pos_strand == 0)
{
recover_UC_Read(overlap_read, &R_INF, overlap_list->list[j].y_id);
@@ -817,9 +1145,11 @@ void ha_print_ovlp_stat(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long
void fill_chain(Fake_Cigar* chain, char* x_string, char* y_string, long long xBeg, long long yBeg,
long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)
{
long long i, xOffset, yOffset, xRegionLen, yRegionLen, /**bandLen,**/ maxXpos, maxYpos, mapScore, zdroped;
/**
long long i, xOffset, yOffset, xRegionLen, yRegionLen, maxXpos, maxYpos, mapGlobalScore, mapExtentScore, zdroped;
long long xBuoundaryScore, yBuoundaryScore;
///float band_rate = 0.08;
int endbouns;
int endbouns,mode;
if(chain->length <= 0) return;
kvec_t(uint8_t) x_num;
@@ -848,7 +1178,9 @@ long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)
///text is x, query is y
afine_gap_alignment(x_string, x_num.a, xRegionLen, y_string, y_num.a, yRegionLen,
c2n, BACKWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW,
/**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapScore, &zdroped);
BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore,
&mapExtentScore, &xBuoundaryScore, &yBuoundaryScore, &zdroped);
// fprintf(stderr, "* xOffset: %lld, yOffset: %lld, xRegionLen: %lld, yRegionLen: %lld, bandLen: %lld, maxXpos: %lld, maxYpos: %lld, zdroped: %lld\n",
// xOffset, yOffset, xRegionLen, yRegionLen, BAND_KSW, maxXpos, maxYpos, zdroped);
}
@@ -889,7 +1221,8 @@ long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)
///text is x, query is y
afine_gap_alignment(x_string+xOffset, x_num.a, xRegionLen, y_string+yOffset, y_num.a, yRegionLen,
c2n, FORWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW,
/**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapScore, &zdroped);
BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore,
&mapExtentScore, &xBuoundaryScore, &yBuoundaryScore, &zdroped);
// fprintf(stderr, "# xOffset: %lld, yOffset: %lld, xRegionLen: %lld, yRegionLen: %lld, bandLen: %lld, maxXpos: %lld, maxYpos: %lld, zdroped: %lld\n",
// xOffset, yOffset, xRegionLen, yRegionLen, BAND_KSW, maxXpos, maxYpos, zdroped);
}
@@ -897,6 +1230,7 @@ long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n)
kv_destroy(x_num);
kv_destroy(y_num);
**/
}
void Final_phasing(overlap_region_alloc* overlap_list, Cigar_record_alloc* cigarline,
UC_Read* g_read, UC_Read* overlap_read, uint8_t* c2n)
@@ -950,20 +1284,10 @@ UC_Read* g_read, UC_Read* overlap_read, uint8_t* c2n)
static void worker_ov_final(void *data, long i, int tid)
{
ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
uint8_t c2n[256]; // this may be moved to ha_ovec_buf_t, but it should be fast to populate anyway
memset(c2n, 4, 256);
c2n[(uint8_t)'A'] = c2n[(uint8_t)'a'] = 0; c2n[(uint8_t)'C'] = c2n[(uint8_t)'c'] = 1;
c2n[(uint8_t)'G'] = c2n[(uint8_t)'g'] = 2; c2n[(uint8_t)'T'] = c2n[(uint8_t)'t'] = 3; // build the encoding table
//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.001, asm_opt.max_n_chain, 0);
/**
correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon,
&matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1,
&current_cigar, &hap, &second_round, 0, 0);
push_final_overlaps(&(R_INF.paf[i]), &overlap_list);
**/
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);
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);
@@ -979,6 +1303,81 @@ static void worker_ov_final(void *data, long i, int tid)
push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
}
void reset_final_overlaps(overlap_region_alloc *overlap_list)
{
uint64_t i;
for (i = 0; i < overlap_list->length; i++)
{
if (overlap_list->list[i].is_match == 1 || overlap_list->list[i].is_match == 2)
{
overlap_list->list[i].x_pos_s = overlap_list->list[i].x_pos_e = (uint32_t)-1;
overlap_list->list[i].y_pos_s = overlap_list->list[i].y_pos_e = (uint32_t)-1;
overlap_list->list[i].is_match = 0;
}
}
ha_sort_list_by_anchor(overlap_list);
}
void debug_affine_gap_alignment(overlap_region_alloc *overlap_list, UC_Read* g_read, UC_Read* overlap_read)
{
uint64_t i;
kvec_t(uint8_t) x_num;
kvec_t(uint8_t) y_num;
kv_init(x_num);
kv_init(y_num);
for (i = 0; i < overlap_list->length; i++)
{
if (overlap_list->list[i].is_match == 1 && overlap_list->list[i].shared_seed == 1)
{
kv_resize(uint8_t, x_num, (uint64_t)(Get_READ_LENGTH(R_INF, overlap_list->list[i].x_id)));
kv_resize(uint8_t, y_num, (uint64_t)(Get_READ_LENGTH(R_INF, overlap_list->list[i].y_id)));
get_affine_gap_score(&(overlap_list->list[i]), g_read, overlap_read, x_num.a, y_num.a,
overlap_list->list[i].x_pos_e + 1 - overlap_list->list[i].x_pos_s,
overlap_list->list[i].y_pos_e + 1 - overlap_list->list[i].y_pos_s);
}
}
kv_destroy(x_num);
kv_destroy(y_num);
}
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);
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);
ma_hit_sort_tn(R_INF.reverse_paf[i].buffer, R_INF.reverse_paf[i].length);
///update_overlaps(&b->olist, &(R_INF.paf[i]), &b->self_read, &b->ovlp_read, 1, 1);
update_overlaps_chain_width(&b->olist, &(R_INF.paf[i]), &b->self_read, &b->ovlp_read, 1, 1, 0.002);
update_overlaps(&b->olist, &(R_INF.reverse_paf[i]), &b->self_read, &b->ovlp_read, 2, 0);
///recover missing exact overlaps
update_exact_overlaps(&b->olist, &b->self_read, &b->ovlp_read);
///Final_phasing(&overlap_list, &cigarline, &g_read, &overlap_read, c2n);
push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, &b->olist, 1);
push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
///debug_affine_gap_alignment(&b->olist, &b->self_read, &b->ovlp_read);
reset_final_overlaps(&b->olist);
correct_overlap_high_het(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read);
push_final_overlaps_increment(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
}
void Output_PAF()
{
fprintf(stderr, "Writing PAF to disk ...... \n");
@@ -1141,41 +1540,99 @@ long long readNum, long long rescue_threshold, float cluster_threshold)
fprintf(stderr, "[M::%s] took %0.2fs, revise edges #: %lld\n\n", __func__, Get_T()-startTime, revises);
}
void hap_recalculate_peaks(char* output_file_name)
{
destory_read_bin(&R_INF);
destory_ma_hit_t_alloc(R_INF.paf);
destory_ma_hit_t_alloc(R_INF.reverse_paf);
char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
sprintf(gfa_name, "%s.ec", output_file_name);
int hom_cov, het_cov;
// construct hash table for high occurrence k-mers
if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0);
ha_opt_update_cov(&asm_opt, hom_cov);
}
free(R_INF.read_length);
free(R_INF.name_index);
load_All_reads(&R_INF, gfa_name);
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
asm_opt.hom_cov = hom_cov;
asm_opt.het_cov = het_cov;
ha_pt_destroy(ha_idx);
ha_idx = 0;
destory_read_bin(&R_INF);
free(gfa_name);
load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name);
fprintf(stderr, "M::%s has done.\n", __func__);
}
void ha_overlap_final(void)
{
int i, hom_cov;
int i, hom_cov, het_cov;
ha_ovec_buf_t **b;
ha_flt_tab_hp = ha_idx_hp = NULL;
CALLOC(b, asm_opt.thread_num);
for (i = 0; i < asm_opt.thread_num; ++i)
b[i] = ha_ovec_init(1, 1);
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov); // build the index
kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads);
b[i] = ha_ovec_init(asm_opt.flag & HA_F_HIGH_HET, 1);///b[i] = ha_ovec_init(1, 1);
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
if(asm_opt.flag & HA_F_HIGH_HET)
{
kt_for(asm_opt.thread_num, worker_ov_final_high_het, b, R_INF.total_reads);
}
else
{
kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads);
}
ha_pt_destroy(ha_idx);
ha_idx = 0;
for (i = 0; i < asm_opt.thread_num; ++i)
ha_ovec_destroy(b[i]);
free(b);
asm_opt.hom_cov = hom_cov;
asm_opt.het_cov = het_cov;
}
int ha_assemble(void)
{
extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
int r, hom_cov = -1, ovlp_loaded = 0;
if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
ovlp_loaded = 1;
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
if (asm_opt.extract_list) {
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), 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;
}
if (!ovlp_loaded) {
ha_flt_tab = ha_idx = NULL;
if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name), load_ct_index(&ha_ct_table, asm_opt.output_file_name);
// construct hash table for high occurrence k-mers
if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov);
if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL)
{
ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0);
ha_opt_update_cov(&asm_opt, hom_cov);
}
// error correction
assert(asm_opt.number_of_round > 0);
for (r = 0; r < asm_opt.number_of_round; ++r) {
for (r = ha_idx?asm_opt.number_of_round-1:0; r < asm_opt.number_of_round; ++r) {
ha_opt_reset_to_round(&asm_opt, r); // this update asm_opt.roundID and a few other fields
ha_overlap_and_correct(r);
fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(),
@@ -1195,9 +1652,11 @@ int ha_assemble(void)
if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
ha_triobin(&asm_opt);
}
build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf,
R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round,
asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded);
if(ovlp_loaded == 2) ovlp_loaded = 0;
build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf,
R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round,
asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded);
destory_All_reads(&R_INF);
return 0;
}
+5
View File
@@ -8,6 +8,11 @@
#define Get_Cigar_Type(RECORD) (RECORD&3)
#define Get_Cigar_Length(RECORD) (RECORD>>2)
#define RESEED_DP 4
#define RESEED_PEAK_RATE 0.15
#define RESEED_LEN 2000
#define RESEED_HP_RATE 0.9
int ha_assemble(void);
#endif
+255 -44
View File
@@ -1,6 +1,8 @@
#include <zlib.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <limits.h>
#include <sys/time.h>
#include "CommandLines.h"
#include "ketopt.h"
@@ -10,11 +12,24 @@
hifiasm_opt_t asm_opt;
static ko_longopt_t long_options[] = {
{ "version", ko_no_argument, 300 },
{ "dbg-gfa", ko_no_argument, 301 },
{ "write-paf", ko_no_argument, 302 },
{ "write-ec", ko_no_argument, 303 },
{ "skip-triobin", ko_no_argument, 304 },
{ "version", ko_no_argument, 300 },
{ "dbg-gfa", ko_no_argument, 301 },
{ "write-paf", ko_no_argument, 302 },
{ "write-ec", ko_no_argument, 303 },
{ "skip-triobin", ko_no_argument, 304 },
{ "max-od-ec", ko_no_argument, 305 },
{ "max-od-final", ko_no_argument, 306 },
{ "ex-list", ko_required_argument, 307 },
{ "ex-iter", ko_required_argument, 308 },
{ "purge-cov", ko_required_argument, 309 },
{ "pri-range", ko_required_argument, 310 },
{ "high-het", ko_no_argument, 311 },
{ "lowQ", ko_required_argument, 312 },
{ "min-hist-cnt", ko_required_argument, 313 },
{ "h1", ko_required_argument, 314 },
{ "h2", ko_required_argument, 315 },
{ "enzyme", ko_required_argument, 316 },
{ "b-cov", ko_required_argument, 317 },
{ 0, 0, 0 }
};
@@ -29,42 +44,56 @@ void Print_H(hifiasm_opt_t* asm_opt)
{
fprintf(stderr, "Usage: hifiasm [options] <in_1.fq> <in_2.fq> <...>\n");
fprintf(stderr, "Options:\n");
fprintf(stderr, " Input/Output:\n");
fprintf(stderr, " -o STR prefix of output files [%s]\n", asm_opt->output_file_name);
fprintf(stderr, " -i ignore saved read correction and overlaps\n");
fprintf(stderr, " -t INT number of threads [%d]\n", asm_opt->thread_num);
fprintf(stderr, " -z INT length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
fprintf(stderr, " --version show version number\n");
fprintf(stderr, " Overlap/Error correction:\n");
fprintf(stderr, " -k INT k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
fprintf(stderr, " -w INT minimizer window size [%d]\n", asm_opt->mz_win);
fprintf(stderr, " -f INT number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift);
fprintf(stderr, " -D FLOAT drop k-mers occurring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
fprintf(stderr, " -N INT consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
fprintf(stderr, " -r INT round of correction [%d]\n", asm_opt->number_of_round);
fprintf(stderr, " Assembly:\n");
fprintf(stderr, " -o FILE prefix of output files [%s]\n", asm_opt->output_file_name);
fprintf(stderr, " -t INT number of threads [%d]\n", asm_opt->thread_num);
fprintf(stderr, " -r INT round of correction [%d]\n", asm_opt->number_of_round);
fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round);
fprintf(stderr, " -k INT k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
fprintf(stderr, " -w INT minimizer window size [%d]\n", asm_opt->mz_win);
fprintf(stderr, " -f INT number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift);
fprintf(stderr, " -D FLOAT drop k-mers occuring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
fprintf(stderr, " -N INT consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
fprintf(stderr, " -i ignore saved overlaps in *.ovlp* files\n");
fprintf(stderr, " -z INT length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
fprintf(stderr, " -m INT size of popped large bubbles for contig graph [%lld]\n", asm_opt->large_pop_bubble_size);
fprintf(stderr, " -p INT size of popped small bubbles for haplotype-resolved unitig graph [%lld]\n", asm_opt->small_pop_bubble_size);
fprintf(stderr, " -n INT small removed unitig threshold [%d]\n", asm_opt->max_short_tip);
fprintf(stderr, " -x FLOAT max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate);
fprintf(stderr, " -y FLOAT min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate);
fprintf(stderr, " --version show version number\n");
fprintf(stderr, " -h show help information\n");
fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round);
fprintf(stderr, " -m INT pop bubbles of <INT in size in contig graphs [%lld]\n", asm_opt->large_pop_bubble_size);
fprintf(stderr, " -p INT pop bubbles of <INT in size in unitig graphs [%lld]\n", asm_opt->small_pop_bubble_size);
fprintf(stderr, " -n INT remove tip unitigs composed of <=INT reads [%d]\n", asm_opt->max_short_tip);
fprintf(stderr, " -x FLOAT max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate);
fprintf(stderr, " -y FLOAT min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate);
fprintf(stderr, " -u disable post join contigs step which may improve N50\n");
fprintf(stderr, " --lowQ INT\n");
fprintf(stderr, " output contig regions with >=INT%% inconsistency in BED format; 0 to disable [%d]\n", asm_opt->bed_inconsist_rate);
fprintf(stderr, " --b-cov INT\n");
fprintf(stderr, " break contigs at breakpoints with coverage drop at <INT-fold coverage [%d]\n", asm_opt->break_cov);
// fprintf(stderr, " --pri-range INT1[,INT2]\n");
// fprintf(stderr, " keep contigs with coverage in this range in p_ctg.gfa; -1 to disable [auto,inf]\n");
fprintf(stderr, " Trio-partition:\n");
fprintf(stderr, " -1 FILE hap1/paternal k-mer dump generated by \"yak count\" []\n");
fprintf(stderr, " -2 FILE hap2/maternal k-mer dump generated by \"yak count\" []\n");
fprintf(stderr, " -3 FILE list of hap1/paternal read names []\n");
fprintf(stderr, " -4 FILE list of hap2/maternal read names []\n");
fprintf(stderr, " -c INT lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt);
fprintf(stderr, " -d INT upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt);
fprintf(stderr, " -1 FILE hap1/paternal k-mer dump generated by \"yak count\" []\n");
fprintf(stderr, " -2 FILE hap2/maternal k-mer dump generated by \"yak count\" []\n");
fprintf(stderr, " -c INT lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt);
fprintf(stderr, " -d INT upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt);
fprintf(stderr, " -3 FILE list of hap1/paternal read names []\n");
fprintf(stderr, " -4 FILE list of hap2/maternal read names []\n");
fprintf(stderr, " Purge-dups:\n");
fprintf(stderr, " -l INT level of purge-dup. In default, [%d] for non-trio; [%d] for trio (see hifiasm.1 for details)\n",
asm_opt->purge_level_primary, asm_opt->purge_level_trio);
fprintf(stderr, " -s FLOAT similarity threshold for duplicate haplotigs [%g]\n",
asm_opt->purge_simi_rate);
fprintf(stderr, " -O INT min number of overlapped reads for duplicate haplotigs [%d]\n",
asm_opt->purge_overlap_len);
fprintf(stderr, " -l INT purge level. 0: no purging; 1: light; 2: aggressive [0 for trio; 2 for unzip]\n");
fprintf(stderr, " -s FLOAT similarity threshold for duplicate haplotigs [%g]\n",
asm_opt->purge_simi_rate);
fprintf(stderr, " -O INT min number of overlapped reads for duplicate haplotigs [%d]\n",
asm_opt->purge_overlap_len);
fprintf(stderr, " --purge-cov INT\n");
fprintf(stderr, " coverage upper bound of Purge-dups [auto]\n");
fprintf(stderr, " --high-het enable this mode for high heterozygosity sample [experimental, not stable]\n");
fprintf(stderr, " Hi-C-partition [experimental, not stable]:\n");
fprintf(stderr, " --h1 FILEs file names of Hi-C R1 [r1_1.fq,r1_2.fq,...]\n");
fprintf(stderr, " --h2 FILEs file names of Hi-C R2 [r2_1.fq,r2_2.fq,...]\n");
fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
fprintf(stderr, "See `man ./hifiasm.1' for detailed description of these command-line options.\n");
@@ -79,15 +108,21 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->read_file_names = NULL;
asm_opt->output_file_name = (char*)(DEFAULT_OUTPUT);
asm_opt->required_read_name = NULL;
asm_opt->hic_enzymes = NULL;
asm_opt->hic_reads[0] = NULL;
asm_opt->hic_reads[1] = NULL;
asm_opt->thread_num = 1;
asm_opt->k_mer_length = 51;
asm_opt->hic_mer_length = 31;
asm_opt->mz_win = 51;
asm_opt->bf_shift = 37;
asm_opt->high_factor = 5.0f;
asm_opt->high_factor = 5.0;
asm_opt->max_ov_diff_ec = 0.04;
asm_opt->max_ov_diff_final = 0.03;
asm_opt->hom_cov = 20;
asm_opt->het_cov = -1024;
asm_opt->max_n_chain = 100;
asm_opt->k_mer_min_freq = 3;
asm_opt->k_mer_max_freq = 66;
asm_opt->min_hist_kmer_cnt = 5;
asm_opt->load_index_from_disk = 1;
asm_opt->write_index_to_disk = 1;
asm_opt->number_of_round = 3;
@@ -108,15 +143,40 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->purge_level_primary = 2;
asm_opt->purge_level_trio = 0;
asm_opt->purge_simi_rate = 0.75;
asm_opt->purge_simi_rate_hic = 0.85;
asm_opt->purge_overlap_len = 1;
asm_opt->purge_overlap_len_hic = 50;
asm_opt->recover_atg_cov_min = -1024;
asm_opt->recover_atg_cov_max = INT_MAX;
asm_opt->hom_global_coverage = -1;
asm_opt->bed_inconsist_rate = 70;
asm_opt->hic_inconsist_rate = 30;
///asm_opt->bub_mer_length = 3;
asm_opt->bub_mer_length = 1000000;
asm_opt->break_cov = 0;
}
void destory_enzyme(enzyme* f)
{
int i;
if(f != NULL)
{
for (i = 0; i < f->n; i++)
{
free(f->a[i]);
}
free(f->a);
free(f->l);
free(f);
}
}
void destory_opt(hifiasm_opt_t* asm_opt)
{
if(asm_opt->read_file_names != NULL)
{
free(asm_opt->read_file_names);
}
if(asm_opt->read_file_names != NULL) free(asm_opt->read_file_names);
if(asm_opt->hic_enzymes != NULL) destory_enzyme(asm_opt->hic_enzymes);
if(asm_opt->hic_reads[0] != NULL) destory_enzyme(asm_opt->hic_reads[0]);
if(asm_opt->hic_reads[1] != NULL) destory_enzyme(asm_opt->hic_reads[1]);
}
void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round)
@@ -156,6 +216,16 @@ static int check_file(char* name, const char* opt)
return 1;
}
static int check_hic_reads(enzyme* f, const char* opt)
{
int i;
for (i = 0; i < f->n; i++)
{
if(check_file(f->a[i], opt) == 0) return 0;
}
return 1;
}
int check_option(hifiasm_opt_t* asm_opt)
{
if(asm_opt->read_file_names == NULL || asm_opt->num_reads == 0)
@@ -264,6 +334,15 @@ int check_option(hifiasm_opt_t* asm_opt)
return 0;
}
if (asm_opt->max_ov_diff_ec < asm_opt->max_ov_diff_final) {
fprintf(stderr, "[ERROR] max_ov_diff_ec shouldn't be smaller than max_ov_diff_final\n");
return 0;
}
if (asm_opt->max_ov_diff_ec < HA_MIN_OV_DIFF) {
fprintf(stderr, "[ERROR] max_ov_diff_ec shouldn't be smaller than %g\n", HA_MIN_OV_DIFF);
return 0;
}
if(asm_opt->max_short_tip < 0)
{
@@ -283,12 +362,62 @@ int check_option(hifiasm_opt_t* asm_opt)
return 0;
}
if(asm_opt->hom_global_coverage < 0 && asm_opt->hom_global_coverage != -1)
{
fprintf(stderr, "[ERROR] purge duplication coverage threshold should be >= 0 (--purge-cov)\n");
return 0;
}
if(asm_opt->bed_inconsist_rate < 0 || asm_opt->bed_inconsist_rate > 100)
{
fprintf(stderr, "[ERROR] inconsistency rate should be [0, 100] (--lowQ)\n");
return 0;
}
if(asm_opt->fn_bin_yak[0] != NULL && check_file(asm_opt->fn_bin_yak[0], "YAK1") == 0) return 0;
if(asm_opt->fn_bin_yak[1] != NULL && check_file(asm_opt->fn_bin_yak[1], "YAK2") == 0) return 0;
if(asm_opt->fn_bin_list[0] != NULL && check_file(asm_opt->fn_bin_list[0], "LIST1") == 0) return 0;
if(asm_opt->fn_bin_list[1] != NULL && check_file(asm_opt->fn_bin_list[1], "LIST2") == 0) return 0;
if(asm_opt->required_read_name != NULL && check_file(asm_opt->required_read_name, "b") == 0) return 0;
if(asm_opt->hic_reads[0] != NULL && check_hic_reads(asm_opt->hic_reads[0], "HIC1") == 0) return 0;
if(asm_opt->hic_reads[1] != NULL && check_hic_reads(asm_opt->hic_reads[1], "HIC2") == 0) return 0;
if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[1] == NULL)
{
fprintf(stderr, "[ERROR] lack r2 of HiC reads (--h2)\n");
return 0;
}
if(asm_opt->hic_reads[1] != NULL && asm_opt->hic_reads[0] == NULL)
{
fprintf(stderr, "[ERROR] lack r1 of HiC reads (--h1)\n");
return 0;
}
if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[1] != NULL &&
asm_opt->hic_reads[0]->n != asm_opt->hic_reads[1]->n)
{
fprintf(stderr, "[ERROR] wrong r1 and r2 of HiC reads (--h1 && --h2)\n");
return 0;
}
if(asm_opt->hic_enzymes != NULL && asm_opt->hic_enzymes->n == 0)
{
fprintf(stderr, "[ERROR] wrong HiC enzymes (--enzyme)\n");
return 0;
}
if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[0]->n == 0)
{
fprintf(stderr, "[ERROR] wrong r1 of HiC reads (--h1)\n");
return 0;
}
if(asm_opt->hic_reads[1] != NULL && asm_opt->hic_reads[1]->n == 0)
{
fprintf(stderr, "[ERROR] wrong r2 of HiC reads (--h2)\n");
return 0;
}
// fprintf(stderr, "input file num: %d\n", asm_opt->num_reads);
// fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
@@ -339,6 +468,61 @@ void get_queries(int argc, char *argv[], ketopt_t* opt, hifiasm_opt_t* asm_opt)
}
}
void get_hic_enzymes(char *argv, enzyme** x, int check_name)
{
int i, k, pre_i, len = strlen(argv);
(*x) = (enzyme*)calloc(1, sizeof(enzyme));
if(len == 0)
{
(*x)->n = 0; (*x)->l = NULL; (*x)->a = NULL;
return;
}
(*x)->n = 1;
for (i = pre_i = 0; i < len; i++)
{
if(argv[i] == ',')
{
(*x)->n++;
continue;
}
if(check_name)
{
if(argv[i] != 'A' && argv[i] != 'C' && argv[i] != 'G' && argv[i] != 'T' &&
argv[i] != 'a' && argv[i] != 'c' && argv[i] != 'g' && argv[i] != 't' &&
argv[i] != 'N' && argv[i] != 'n')
{
(*x)->n = 0;
(*x)->l = NULL;
(*x)->a = NULL;
return;
}
}
}
(*x)->l = (int*)calloc((*x)->n, sizeof(int));
(*x)->a = (char**)calloc((*x)->n, sizeof(char*));
for (i = pre_i = k = 0; i < len; i++)
{
if(argv[i] == ',')
{
(*x)->l[k] = i - pre_i;
(*x)->a[k] = (char*)malloc(sizeof(char)*((*x)->l[k]+1));
memcpy((*x)->a[k], argv + pre_i, (*x)->l[k]);
(*x)->a[k][(*x)->l[k]] = '\0';
pre_i = i + 1;
k++;
}
}
(*x)->l[k] = i - pre_i;
(*x)->a[k] = (char*)malloc(sizeof(char)*((*x)->l[k]+1));
memcpy((*x)->a[k], argv + pre_i, (*x)->l[k]);
(*x)->a[k][(*x)->l[k]] = '\0';
}
int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
{
@@ -346,7 +530,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
int c;
while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:w:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:1:2:3:4:l:s:O:", long_options)) >= 0) {
while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:w:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:1:2:3:4:l:s:O:eu", long_options)) >= 0) {
if (c == 'h')
{
Print_H(asm_opt);
@@ -381,10 +565,34 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 'p') asm_opt->small_pop_bubble_size = atoll(opt.arg);
else if (c == 'm') asm_opt->large_pop_bubble_size = atoll(opt.arg);
else if (c == 'n') asm_opt->max_short_tip = atoll(opt.arg);
else if (c == 'e') asm_opt->flag |= HA_F_BAN_ASSEMBLY;
else if (c == 'u') asm_opt->flag |= HA_F_BAN_POST_JOIN;
else if (c == 301) asm_opt->flag |= HA_F_VERBOSE_GFA;
else if (c == 302) asm_opt->flag |= HA_F_WRITE_PAF;
else if (c == 303) asm_opt->flag |= HA_F_WRITE_EC;
else if (c == 304) asm_opt->flag |= HA_F_SKIP_TRIOBIN;
else if (c == 305) asm_opt->max_ov_diff_ec = atof(opt.arg);
else if (c == 306) asm_opt->max_ov_diff_final = atof(opt.arg);
else if (c == 307) asm_opt->extract_list = opt.arg;
else if (c == 308) asm_opt->extract_iter = atoi(opt.arg);
else if (c == 309) asm_opt->hom_global_coverage = atoi(opt.arg);
else if (c == 310)
{
char* s = NULL;
asm_opt->recover_atg_cov_min = strtol(opt.arg, &s, 10);
if (*s == ',') asm_opt->recover_atg_cov_max = strtol(s + 1, &s, 10);
if(asm_opt->recover_atg_cov_min == -1 || asm_opt->recover_atg_cov_max == -1)
{
asm_opt->recover_atg_cov_min = asm_opt->recover_atg_cov_max = -1;
}
}
else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET;
else if (c == 312) asm_opt->bed_inconsist_rate = atoi(opt.arg);
else if (c == 313) asm_opt->min_hist_kmer_cnt = atoi(opt.arg);
else if (c == 314) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[0]), 0);
else if (c == 315) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[1]), 0);
else if (c == 316) get_hic_enzymes(opt.arg, &(asm_opt->hic_enzymes), 1);
else if (c == 317) asm_opt->break_cov = atoi(opt.arg);
else if (c == 'l')
{ ///0: disable purge_dup; 1: purge containment; 2: purge overlap
asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);
@@ -403,6 +611,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
}
}
if (argc == opt.ind)
{
Print_H(asm_opt);
@@ -411,5 +620,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
get_queries(argc, argv, &opt, asm_opt);
return check_option(asm_opt);
}
+43 -11
View File
@@ -3,18 +3,28 @@
#include <pthread.h>
#define HA_VERSION "0.5-dirty-r247"
#define HA_VERSION "0.14-r309"
#define VERBOSE 0
#define HA_F_NO_HPC 0x1
#define HA_F_NO_KMER_FLT 0x2
#define HA_F_VERBOSE_GFA 0x4
#define HA_F_WRITE_EC 0x8
#define HA_F_WRITE_PAF 0x10
#define HA_F_SKIP_TRIOBIN 0x20
#define HA_F_NO_HPC 0x1
#define HA_F_NO_KMER_FLT 0x2
#define HA_F_VERBOSE_GFA 0x4
#define HA_F_WRITE_EC 0x8
#define HA_F_WRITE_PAF 0x10
#define HA_F_SKIP_TRIOBIN 0x20
#define HA_F_PURGE_CONTAIN 0x40
#define HA_F_PURGE_JOIN 0x80
#define HA_F_PURGE_JOIN 0x80
#define HA_F_BAN_POST_JOIN 0x100
#define HA_F_BAN_ASSEMBLY 0x200
#define HA_F_HIGH_HET 0x400
#define HA_MIN_OV_DIFF 0.02 // min sequence divergence in an overlap
typedef struct{
int *l, n;
char **a;
}enzyme;
typedef struct {
int flag;
@@ -24,15 +34,24 @@ typedef struct {
char* required_read_name;
char *fn_bin_yak[2];
char *fn_bin_list[2];
char *extract_list;
enzyme *hic_reads[2];
enzyme *hic_enzymes;
int extract_iter;
int thread_num;
int k_mer_length;
int hic_mer_length;
int bub_mer_length;
int mz_win;
int bf_shift;
float high_factor; // coverage cutoff set to high_factor*hom_cov
double high_factor; // coverage cutoff set to high_factor*hom_cov
double max_ov_diff_ec;
double max_ov_diff_final;
int hom_cov;
int het_cov;
int break_cov;
int max_n_chain; // fall-back max number of chains to consider
int k_mer_min_freq;
int k_mer_max_freq;
int min_hist_kmer_cnt;
int load_index_from_disk;
int write_index_to_disk;
int number_of_round;
@@ -49,11 +68,18 @@ typedef struct {
int purge_level_primary;
int purge_level_trio;
int purge_overlap_len;
int purge_overlap_len_hic;
int recover_atg_cov_min;
int recover_atg_cov_max;
int hom_global_coverage;
int bed_inconsist_rate;
int hic_inconsist_rate;
float max_hang_rate;
float min_drop_rate;
float max_drop_rate;
float purge_simi_rate;
float purge_simi_rate_hic;
long long small_pop_bubble_size;
long long large_pop_bubble_size;
@@ -62,6 +88,7 @@ typedef struct {
long long num_recorrected_bases;
long long mem_buf;
long long coverage;
} hifiasm_opt_t;
extern hifiasm_opt_t asm_opt;
@@ -78,4 +105,9 @@ static inline int ha_opt_triobin(const hifiasm_opt_t *opt)
return ((opt->fn_bin_yak[0] && opt->fn_bin_yak[1]) || (opt->fn_bin_list[0] && opt->fn_bin_list[1]));
}
static inline int ha_opt_hic(const hifiasm_opt_t *opt)
{
return ((opt->hic_reads[0] && opt->hic_reads[1]));
}
#endif
+1014 -50
View File
File diff suppressed because it is too large Load Diff
+30 -6
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
@@ -16,11 +17,12 @@
#define INSERTION 2
#define DELETION 3
#define WINDOW_MAX_SIZE (WINDOW + (int)(1.0 / HA_MIN_OV_DIFF) + 3) // TODO: why 1/max_ov_diff?
///#define FLAG_THRE 0
#define MAX(x, y) ((x >= y)?(x):(y))
#define MIN(x, y) ((x <= y)?(x):(y))
#define MAX(x, y) (((x) >= (y))?(x):(y))
#define MIN(x, y) (((x) <= (y))?(x):(y))
#define DIFF(x, y) ((MAX((x), (y))) - (MIN((x), (y))))
#define OVERLAP(x_start, x_end, y_start, y_end) (MIN(x_end, y_end) - MAX(x_start, y_start) + 1)
///#define OVERLAP(x_start, x_end, y_start, y_end) MIN(x_end, y_end) - MAX(x_start, y_start) + 1
@@ -1160,10 +1162,32 @@ void init_Cigar_record_alloc(Cigar_record_alloc* x);
void resize_Cigar_record_alloc(Cigar_record_alloc* x, long long new_size);
void destory_Cigar_record_alloc(Cigar_record_alloc* x);
void afine_gap_alignment(const char *tseq, uint8_t* tnum, const int tl,
const char *qseq, uint8_t* qnum, const int ql, const uint8_t *c2n, const int strand,
void afine_gap_alignment(const char *qseq, uint8_t* qnum, const int ql,
const char *tseq, uint8_t* tnum, const int tl, const uint8_t *c2n, const int strand,
int sc_mch, int sc_mis, int gapo, int gape, int bandLen, int zdrop, int end_bonus,
long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped);
long long* max_q_pos, long long* max_t_pos, long long* global_score,
long long* extention_score, long long* q_boundary_score, long long* q_boundary_t_coordinate,
long long* t_boundary_score, long long* t_boundary_q_coordinate,
long long* droped, int mode);
void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_INF,
UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read);
long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* overlap_read, uint8_t* x_num,
uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen);
int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u8_warp* k_flag, float hp_rate, int rlen, FILE* fp);
inline int if_exact_match(char* x, long long xLen, char* y, long long yLen, long long xBeg, long long xEnd, long long yBeg, long long yEnd)
{
long long overlapLen = xEnd - xBeg + 1;
if(yEnd - yBeg + 1 == overlapLen)
{
if(memcmp(x+xBeg, y+yBeg, overlapLen)==0) return 1;
}
return 0;
}
#define FORWARD_KSW 0
#define BACKWARD_KSW 1
@@ -1172,5 +1196,5 @@ long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped)
#define GAP_OPEN_KSW 4
#define GAP_EXT_KSW 2
#define Z_DROP_KSW 400
#define BAND_KSW 50
#define BAND_KSW 500
#endif
+125 -5
View File
@@ -305,6 +305,9 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp)
}
}
long long get_chainLen(long long x_beg, long long x_end, long long xLen,
long long y_beg, long long y_end, long long yLen)
{
@@ -337,6 +340,50 @@ long long get_chainLen(long long x_beg, long long x_end, long long xLen,
return x_end - x_beg + 1;
}
void debug_chain_single_site(k_mer_hit* a, long long a_n, Chain_Data* dp, int x_readLen, int y_readLen, int s_index)
{
long long j, current_j = s_index;
long long selfLen = 0, indels = 0;
long long distance_self_pos, distance_pos, distance_gap;
j = s_index;
while (j >= 0)
{
current_j = j;
j = dp->pre[j];
if(j != -1)
{
distance_self_pos = a[current_j].self_offset - a[j].self_offset;
distance_pos = a[current_j].offset - a[j].offset;
distance_gap = distance_pos > distance_self_pos? distance_pos - distance_self_pos : distance_self_pos - distance_pos;
indels += distance_gap;
selfLen += distance_self_pos;
}
fprintf(stderr, "j: %lld, score: %lld, occ: %d, pre_j: %lld\n",
current_j, (long long)dp->score[current_j], dp->occ[current_j], j);
}
fprintf(stderr, "s_self_offset: %u, s_offset: %u, e_self_offset: %u, e_offset: %u, ovlp length: %lld, x_readLen: %d, y_readLen: %d\n",
a[s_index].self_offset, a[s_index].offset, a[current_j].self_offset, a[current_j].offset,
get_chainLen(a[s_index].self_offset, a[current_j].self_offset, x_readLen,
a[s_index].offset, a[current_j].offset, y_readLen), x_readLen, y_readLen);
if(indels != dp->indels[s_index])
{
fprintf(stderr, "indels: %lld, dp->indels[i]: %ld\n", indels, (long)dp->indels[s_index]);
}
if(selfLen != dp->self_length[s_index])
{
fprintf(stderr, "selfLen: %lld, dp->self_length[i]: %ld\n", selfLen, (long)dp->self_length[s_index]);
}
fprintf(stderr,"\n");
}
int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc, double bw_thres)
{
int32_t i, tot_indel = 0, tot_len = 0;
@@ -348,7 +395,7 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc
if (i < n_a) return -1;
bw_pen = 1.0 / bw_thres;
dp->score[0] = a[0].good? min_sc : min_sc>>1;
dp->pre[0] = -1, dp->indels[0] = 0, dp->self_length[0] = 0;
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;
int32_t dx = (int32_t)a[i].offset - (int32_t)a[i-1].offset;
@@ -368,6 +415,7 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc
dp->pre[i] = i - 1;
dp->indels[i] = tot_indel;
dp->self_length[i] = tot_len;
dp->occ[i] = i + 1;
}
if (i < n_a) return -1;
return n_a;
@@ -391,6 +439,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
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;
@@ -435,7 +484,9 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
///min distance
distance_min = distance_pos < distance_self_pos? distance_pos:distance_self_pos;
score = distance_min < min_score? distance_min : min_score;
if (!a[j].good) score >>= 1;
///need to be fixed in r305
///if (!a[j].good) score = (score >> 1) + (score & 1);
if (!a[j].good) score >>= 1;
gap_rate = (double)((double)(total_indels)/(double)(total_self_length));
///if the gap rate > 0.06, score will be negative
@@ -444,7 +495,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
score += dp->score[j];
///find a new max score
if (score > max_score) {
if (score > max_score) {///must use > instead of >=
max_score = score;
max_j = j;
max_indels = total_indels;
@@ -466,10 +517,12 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
dp->pre[i] = max_j;
dp->indels[i] = max_indels;
dp->self_length[i] = max_self_length;
dp->occ[i] = 1;
if(max_j != -1) dp->occ[i] = dp->occ[max_j] + 1;
}
///debug_chain(a, a_n, dp);
skip_dp:
max_score = -1;
@@ -562,7 +615,7 @@ skip_dp:
}
}
void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list,
void calculate_overlap_region_by_chaining_back(Candidates_list* candidates, overlap_region_alloc* overlap_list,
uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end)
{
overlap_region tmp_region;
@@ -627,6 +680,71 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r
destory_fake_cigar(&(tmp_region.f_cigar));
}
void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx,
uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end, overlap_region* f_cigar)
{
long long i = 0;
uint64_t current_ID;
uint64_t current_stand;
if (candidates->length == 0)
{
return;
}
long long sub_region_beg;
long long sub_region_end;
clear_fake_cigar(&((*f_cigar).f_cigar));
i = 0;
while (i < candidates->length)
{
chain_idx->a.n = 0;
current_ID = candidates->list[i].readID;
current_stand = candidates->list[i].strand;
///reference read
(*f_cigar).x_id = readID;
(*f_cigar).x_pos_strand = current_stand;
///query read
(*f_cigar).y_id = current_ID;
///here the strand of query is always 0
(*f_cigar).y_pos_strand = 0;
sub_region_beg = i;
sub_region_end = i;
i++;
while (i < candidates->length
&&
current_ID == candidates->list[i].readID
&&
current_stand == candidates->list[i].strand)
{
sub_region_end = i;
i++;
}
if ((*f_cigar).x_id == (*f_cigar).y_id)
{
continue;
}
chain_DP(candidates->list + sub_region_beg,
sub_region_end - sub_region_beg + 1, &(candidates->chainDP), f_cigar, band_width_threshold,
25, Get_READ_LENGTH((*R_INF), (*f_cigar).x_id), Get_READ_LENGTH((*R_INF), (*f_cigar).y_id));
///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1)
if ((*f_cigar).x_id != (*f_cigar).y_id)
{
append_inexact_overlap_region_alloc(overlap_list, f_cigar, R_INF, add_beg_end);
}
}
}
void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error,
int extra_begin, int extra_end, int error_threshold)
{
@@ -704,6 +822,7 @@ void destory_Chain_Data(Chain_Data* x)
free(x->pre);
free(x->indels);
free(x->self_length);
free(x->occ);
free(x->tmp);
}
@@ -716,6 +835,7 @@ void resize_Chain_Data(Chain_Data* x, long long size)
REALLOC(x->pre, x->size);
REALLOC(x->indels, x->size);
REALLOC(x->self_length, x->size);
REALLOC(x->occ, x->size);
REALLOC(x->tmp, x->size);
}
}
+5 -11
View File
@@ -6,22 +6,15 @@
#define MAX_SUFFIX_BITS 64
#define MODE_VALUE 101
///#define WINDOW 350
///#define THRESHOLD 14
#define WINDOW 375
//#define WINDOW_BOUNDARY 150
#define WINDOW_BOUNDARY 375
///for one side, the first or last WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY bases should not be corrected
#define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25
#define THRESHOLD 15
#define THRESHOLD_RATE 0.04
#define TAIL_LENGTH int(1/THRESHOLD_RATE)
///#define OVERLAP_THRESHOLD 0.9
#define OVERLAP_THRESHOLD_FILTER 0.9
#define WINDOW_MAX_SIZE WINDOW + TAIL_LENGTH + 3
#define HIGH_HET_OVERLAP_THRESHOLD_FILTER 0.3
#define HIGH_HET_ERROR_RATE 0.08
#define THRESHOLD_MAX_SIZE 31
#define FINAL_OVERLAP_ERROR_RATE 0.03
#define GROUP_SIZE 4
///the max cigar likes 10M10D10M10D10M
@@ -126,6 +119,7 @@ typedef struct {
int64_t *pre;
int32_t *indels;
int32_t *self_length;
int32_t *occ;
int64_t *tmp; // MUST BE 64-bit integer
int64_t length;
int64_t size;
@@ -151,8 +145,8 @@ int extra_begin, int extra_end, int error_threshold);
void overlap_region_sort_y_id(overlap_region *a, long long n);
void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list,
uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end);
void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx,
uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end, overlap_region* f_cigar);
void init_fake_cigar(Fake_Cigar* x);
void destory_fake_cigar(Fake_Cigar* x);
+12 -3
View File
@@ -1,10 +1,12 @@
CXX= g++
CC= gcc
CXXFLAGS= -g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall
CFLAGS= $(CXXFLAGS)
CPPFLAGS=
INCLUDES=
OBJS= CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \
htab.o hist.o sketch.o anchor.o sys.o
htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o
EXE= hifiasm
LIBS= -lz -lpthread -lm
@@ -13,12 +15,15 @@ ifneq ($(asan),)
LIBS+=-fsanitize=address
endif
.SUFFIXES:.cpp .o
.SUFFIXES:.cpp .c .o
.PHONY:all clean depend
.cpp.o:
$(CXX) -c $(CXXFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@
.c.o:
$(CC) -c $(CFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@
all:$(EXE)
$(EXE):$(OBJS) main.o
@@ -38,6 +43,7 @@ Assembly.o: kthread.h
CommandLines.o: CommandLines.h ketopt.h
Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h
Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h
Correct.o: ksw2.h
Hash_Table.o: Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h
Hash_Table.o: CommandLines.h ksort.h
Levenshtein_distance.o: Levenshtein_distance.h
@@ -51,10 +57,12 @@ Process_Read.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
Purge_Dups.o: ksort.h Purge_Dups.h kvec.h kdq.h Overlaps.h Hash_Table.h
Purge_Dups.o: htab.h Process_Read.h CommandLines.h Correct.h
Purge_Dups.o: Levenshtein_distance.h POA.h kthread.h
Trio.o: khashl.h kthread.h Process_Read.h Overlaps.h kvec.h kdq.h
Trio.o: khashl.h kthread.h kseq.h Process_Read.h Overlaps.h kvec.h kdq.h
Trio.o: CommandLines.h htab.h
anchor.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
anchor.o: ksort.h Hash_Table.h
extract.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h khashl.h
extract.o: kseq.h
hist.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
htab.o: kthread.h khashl.h kseq.h ksort.h htab.h Process_Read.h Overlaps.h
htab.o: kvec.h kdq.h CommandLines.h
@@ -63,3 +71,4 @@ main.o: CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h Assembly.h
main.o: Levenshtein_distance.h htab.h
sketch.o: kvec.h htab.h Process_Read.h Overlaps.h kdq.h CommandLines.h
sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
hic.o: hic.h
+4838 -4952
View File
File diff suppressed because it is too large Load Diff
+89 -37
View File
@@ -21,11 +21,16 @@
#define PRIMARY_LABLE 0
#define ALTER_LABLE 1
#define HAP_LABLE 2
#define FAKE_LABLE 4
#define TRIO_THRES 0.9
#define DOUBLE_CHECK_THRES 0.1
#define FINAL_DOUBLE_CHECK_THRES 0.2
#define CHIMERIC_TRIM_THRES 4
// #define PRIMARY_LABLE 1
// #define ALTER_LABLE 2
// #define HAP_LABLE 4
#define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32))
#define Get_qs(RECORD) ((uint32_t)((RECORD).qns))
#define Get_qe(RECORD) ((RECORD).qe)
@@ -102,6 +107,11 @@ typedef struct {
uint8_t no_l_indel;
} asg_arc_t;
typedef struct {
size_t n, m;
asg_arc_t* a;
} kv_asg_arc_t;
typedef struct {
uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero
@@ -151,6 +161,7 @@ typedef struct { size_t n, m; ma_utg_t *a; } ma_utg_v;
typedef struct {
ma_utg_v u;
asg_t *g;
kvec_t(uint64_t) occ;
} ma_ug_t;
typedef struct {
@@ -390,7 +401,8 @@ typedef struct {
}kvec_asg_arc_t_warp;
void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend);
int asg_arc_del_multi(asg_t *g);
int asg_arc_del_asymm(asg_t *g);
typedef struct {
uint32_t q_pos;
@@ -438,11 +450,8 @@ long long max_hang_length, long long clean_round, long long gap_fuzz,
float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name,
long long bubble_dist, int read_graph, int write);
void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, int id, char* command);
void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int id, char* command);
void collect_abnormal_edges(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum);
void add_overlaps(ma_hit_t_alloc* source_paf, ma_hit_t_alloc* dest_paf, uint64_t* source_index, long long listLen);
void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen);
void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf,
@@ -451,18 +460,6 @@ uint64_t* source_index, long long listLen);
#define EvaluateLen(U, id) ((U).a[(id)].start)
#define IsMerge(U, id) ((U).a[(id)].end)
#define kv_reuse(v, rn, rm, r) ((v).n = (rn), (v).m = (rm), (v).a = (r))
#define long_tip(U, id, threshold) ((EvaluateLen((U), (id))>=(threshold))&&(!((U).a[(id)].circ)))
///there are threee cases:
///1. if this untig is too long (>maxShortUntig), it must be not short untig/must be a long untig
///2. if this untig is long (>minLongUntig && EvaluateLen(ug->u, av[i].v>>1) > (EvaluateLen(ug->u, v>>1)*l_untig_rate)), it might be a long tip
#define check_long_tip(U, id, minLongUntig, maxShortUntig, ShortUntigRate, mainLen) \
((!((U).a[(id)].circ)) \
&& \
((EvaluateLen((U), (id)) > (maxShortUntig))\
||\
((long_tip((U), (id), (minLongUntig)))\
&&\
(EvaluateLen((U), (id)) > (ShortUntigRate)*(mainLen)))))
#define Get_vis(visit, v, d) (((visit)[(v)>>1])&(((((v)<<(d))&1)+1)))
#define Set_vis(visit, v, d) (((visit)[(v)>>1])|=(((((v)<<(d))&1)+1)))
@@ -477,14 +474,10 @@ typedef struct {
void init_R_to_U(R_to_U* x, uint64_t len);
void destory_R_to_U(R_to_U* x);
void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig);
void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8_t* flag);
void get_R_to_U(R_to_U* x, uint32_t rID, uint32_t* uID, uint32_t* is_Unitig);
void transfor_R_to_U(R_to_U* x);
void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, int require_equal_nv, int test_tangle);
void clean_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold,
R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen,
uint32_t miniBiGraph, float chimeric_rate, int is_final_clean);
int asg_pop_bubble_primary(asg_t *g, int max_dist);
long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *g, long long circleLen, int is_drop);
@@ -739,13 +732,12 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
stops_threshold, b_0) == LOOP)
{
return UNAVAILABLE;
}
}
if(get_unitig(nsg, ug, v_1, &vEnd, &ELen_1, &tmp, &max_stop_nodeLen, &max_stop_baseLen,
stops_threshold, b_1) == LOOP)
{
return UNAVAILABLE;
}
if(ELen_0<=min_edge_length || ELen_1<=min_edge_length) return UNAVAILABLE;
rIdContig b_max, b_min;
@@ -766,7 +758,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
uint32_t max_count = 0, min_count = 0;
ma_utg_t *node_min = NULL, *node_max = NULL;
if(ug != NULL)
{
/*****************************label all unitigs****************************************/
@@ -777,12 +768,11 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
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);
set_R_to_U(ruIndex, qn, (b_max.b_0->b.a[b_max.untigI]>>1), 1, &(read_sg->seq[qn].c));
}
}
/*****************************label all unitigs****************************************/
///each unitig
for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
{
@@ -816,7 +806,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
}
}
}
/*****************************label all unitigs****************************************/
for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
{
@@ -829,7 +818,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
}
}
/*****************************label all unitigs****************************************/
}
else
{
@@ -837,7 +825,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
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);
set_R_to_U(ruIndex, qn, 1, 1, &(read_sg->seq[qn].c));
}
/*****************************label all reads****************************************/
@@ -892,8 +880,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
return NON_PLOID;
}
inline uint32_t check_different_haps_naive(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg,
uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1,
R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
@@ -1043,10 +1029,8 @@ uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp*
void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t is_bubble_check,
uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges);
void deduplicate(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig,
long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t resolve_tangle);
void all_to_all_deduplicate(ma_ug_t* ug, uint8_t postive_flag, float drop_rate,
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
void 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);
void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges);
@@ -1057,9 +1041,77 @@ uint32_t is_bubble_check, uint32_t is_primary_check);
uint32_t get_edge_from_source(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t query, uint32_t target, asg_arc_t* t);
uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, int max_dist, buf_t *b,
uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop);
uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes);
int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio);
typedef struct{
double weight;
uint32_t uID:31, del:1;
uint64_t dis;
uint64_t occ;
///uint64_t occ:63, scaff:1;
///uint32_t enzyme;
} hc_edge;
typedef struct{
kvec_t(hc_edge) e;
kvec_t(hc_edge) f;//forbiden
} hc_linkeage;
typedef struct{
uint64_t beg, end;
}bed_interval;
typedef struct{
size_t n, m;
bed_interval* a;
}bed_in;
typedef struct{
kvec_t(hc_linkeage) a;
kvec_t(uint64_t) enzymes;
kvec_t(bed_in) bed;
uint32_t* u_idx;
uint64_t r_num;
} hc_links;
typedef struct{
///kvec_t(hc_edge) a;
size_t n, m;
hc_edge *a;
}hc_edge_warp;
void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num);
void destory_hc_links(hc_links* link);
void clean_primary_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold,
R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen,
uint32_t miniBiGraph, float chimeric_rate, int is_final_clean, int just_bubble_pop,
float drop_ratio, hc_links* link);
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,
R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link);
void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg);
ma_ug_t* copy_untig_graph(ma_ug_t *src);
ma_ug_t* output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name,
uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist,
long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex,
float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench);
asg_t* copy_read_graph(asg_t *src);
ma_ug_t *ma_ug_gen(asg_t *g);
void ma_ug_destroy(ma_ug_t *ug);
inline int inter_interval(int a_s, int a_e, int b_s, int b_e, int* i_s, int* i_e)
{
if(a_s > b_e || b_s > a_e) return 0;
if(i_s) (*i_s) = a_s >= b_s? a_s : b_s; ///MAX(a_s, b_s);
if(i_e) (*i_e) = a_e <= b_e? a_e : b_e; ///MIN(a_e, b_e);
return 1;
}
#define JUNK_COV 5
#define DISCARD_RATE 0.8
+189 -87
View File
@@ -46,6 +46,7 @@ void destory_All_reads(All_reads* r)
if (r->read_sperate[i]) free(r->read_sperate[i]);
if (r->paf && r->paf[i].buffer) free(r->paf[i].buffer);
if (r->reverse_paf && r->reverse_paf[i].buffer) free(r->reverse_paf[i].buffer);
///if (r->pb_regions) kv_destroy(r->pb_regions[i].a);
}
free(r->paf);
free(r->reverse_paf);
@@ -55,6 +56,7 @@ void destory_All_reads(All_reads* r)
free(r->name_index);
free(r->read_length);
free(r->trio_flag);
///if (r->pb_regions) free(r->pb_regions);
}
void write_All_reads(All_reads* r, char* read_file_name)
@@ -98,6 +100,9 @@ void write_All_reads(All_reads* r, char* read_file_name)
fwrite(r->name, sizeof(char), r->total_name_length, fp);
fwrite(r->name_index, sizeof(uint64_t), r->name_index_size, fp);
fwrite(r->trio_flag, sizeof(uint8_t), r->total_reads, fp);
fwrite(&(asm_opt.hom_cov), sizeof(asm_opt.hom_cov), 1, fp);
fwrite(&(asm_opt.het_cov), sizeof(asm_opt.het_cov), 1, fp);
free(index_name);
fflush(fp);
fclose(fp);
@@ -173,6 +178,8 @@ int load_All_reads(All_reads* r, char* read_file_name)
/****************************may have bugs********************************/
r->trio_flag = (uint8_t*)malloc(sizeof(uint8_t)*r->total_reads);
f_flag += fread(r->trio_flag, sizeof(uint8_t), r->total_reads, fp);
f_flag += fread(&(asm_opt.hom_cov), sizeof(asm_opt.hom_cov), 1, fp);
f_flag += fread(&(asm_opt.het_cov), sizeof(asm_opt.het_cov), 1, fp);
/****************************may have bugs********************************/
r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
@@ -187,6 +194,7 @@ int load_All_reads(All_reads* r, char* read_file_name)
r->second_round_cigar[i].lost_base_length = r->cigars[i].lost_base_length = 0;
r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL;
}
///r->pb_regions = NULL;
free(index_name);
fclose(fp);
@@ -195,6 +203,35 @@ int load_All_reads(All_reads* r, char* read_file_name)
return 1;
}
int destory_read_bin(All_reads* r)
{
uint64_t i = 0;
for (i = 0; i < r->total_reads; i++)
{
if (r->N_site[i]) free(r->N_site[i]);
if (r->read_sperate[i]) free(r->read_sperate[i]);
if (r->cigars[i].record) free(r->cigars[i].record);
if (r->cigars[i].lost_base) free(r->cigars[i].lost_base);
if (r->second_round_cigar[i].record) free(r->second_round_cigar[i].record);
if (r->second_round_cigar[i].lost_base) free(r->second_round_cigar[i].lost_base);
}
free(r->N_site);
free(r->read_length);
free(r->read_size);
free(r->read_sperate);
free(r->name);
free(r->name_index);
free(r->trio_flag);
free(r->cigars);
free(r->second_round_cigar);
return 1;
}
void ha_insert_read_len(All_reads *r, int read_len, int name_len)
{
r->total_reads++;
@@ -229,6 +266,8 @@ void malloc_All_reads(All_reads* r)
r->second_round_cigar = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
///r->pb_regions = (kvec_t_u64_warp*)malloc(r->total_reads*sizeof(kvec_t_u64_warp));
for (i = 0; i < (long long)r->total_reads; i++)
{
r->second_round_cigar[i].size = r->cigars[i].size = 0;
@@ -240,6 +279,7 @@ void malloc_All_reads(All_reads* r)
r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL;
init_ma_hit_t_alloc(&(r->paf[i]));
init_ma_hit_t_alloc(&(r->reverse_paf[i]));
///kv_init(r->pb_regions[i].a);
}
r->name = (char*)malloc(sizeof(char)*r->total_name_length);
@@ -298,93 +338,6 @@ void init_UC_Read(UC_Read* r)
}
}
void recover_UC_Read_sub_region_begin_end(char* r, long long start_pos, long long length, uint8_t strand,
All_reads* R_INF, long long ID, int extra_begin, int extra_end)
{
long long readLen = Get_READ_LENGTH((*R_INF), ID);
uint8_t* src = Get_READ((*R_INF), ID);
long long i;
long long copyLen;
long long end_pos = start_pos + length - 1;
if (strand == 0)
{
i = start_pos;
copyLen = 0;
long long initLen = start_pos % 4;
if (initLen != 0)
{
memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen);
copyLen = copyLen + 4 - initLen;
i = i + copyLen;
}
while (copyLen < length)
{
memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4);
copyLen = copyLen + 4;
i = i + 4;
}
if (R_INF->N_site[ID])
{
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
{
if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos)
{
r[R_INF->N_site[ID][i] - start_pos] = 'N';
}
else if((long long)R_INF->N_site[ID][i] > end_pos)
{
break;
}
}
}
}
else
{
start_pos = readLen - start_pos - 1;
end_pos = readLen - end_pos - 1;
///start_pos > end_pos
i = start_pos;
copyLen = 0;
long long initLen = (start_pos + 1) % 4;
if (initLen != 0)
{
memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen);
copyLen = copyLen + initLen;
i = i - initLen;
}
while (copyLen < length)
{
memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4);
copyLen = copyLen + 4;
i = i - 4;
}
if (R_INF->N_site[ID])
{
long long offset = readLen - start_pos - 1;
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
{
if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos)
{
r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N';
}
else if((long long)R_INF->N_site[ID][i] > start_pos)
{
break;
}
}
}
}
}
void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID)
{
@@ -475,6 +428,104 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length,
}
void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID)
{
i_r->length = length;i_r->RID = ID;
if (i_r->length + 8 > i_r->size)
{
i_r->size = i_r->length + 4;
i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size));
}
char* r = i_r->seq;
long long readLen = Get_READ_LENGTH((*R_INF), ID);
uint8_t* src = Get_READ((*R_INF), ID);
long long i;
long long copyLen;
long long end_pos = start_pos + length - 1;
if (strand == 0)
{
i = start_pos;
copyLen = 0;
long long initLen = start_pos % 4;
if (initLen != 0)
{
memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen);
copyLen = copyLen + 4 - initLen;
i = i + copyLen;
}
while (copyLen < length)
{
memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4);
copyLen = copyLen + 4;
i = i + 4;
}
if (R_INF->N_site[ID])
{
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
{
if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos)
{
r[R_INF->N_site[ID][i] - start_pos] = 'N';
}
else if((long long)R_INF->N_site[ID][i] > end_pos)
{
break;
}
}
}
}
else
{
start_pos = readLen - start_pos - 1;
end_pos = readLen - end_pos - 1;
///start_pos > end_pos
i = start_pos;
copyLen = 0;
long long initLen = (start_pos + 1) % 4;
if (initLen != 0)
{
memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen);
copyLen = copyLen + initLen;
i = i - initLen;
}
while (copyLen < length)
{
memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4);
copyLen = copyLen + 4;
i = i - 4;
}
if (R_INF->N_site[ID])
{
long long offset = readLen - start_pos - 1;
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
{
if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos)
{
r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N';
}
else if((long long)R_INF->N_site[ID][i] > start_pos)
{
break;
}
}
}
}
}
void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID)
{
r->length = Get_READ_LENGTH((*R_INF), ID);
@@ -633,3 +684,54 @@ void reverse_complement(char* pattern, uint64_t length)
pattern[end] = RC_CHAR(pattern[end]);
}
}
void init_Debug_reads(Debug_reads* x, const char* file)
{
int nameLen, i, bufLen = 1000;
if((uint64_t)(bufLen) < strlen(file) + 50) bufLen = strlen(file) + 50;
char* Name_Buffer = (char*)malloc(sizeof(char)*bufLen);
fprintf(stderr, "Queried debugging reads at: %s\n", file);
x->fp = fopen(file,"r");
x->query_num = 0;
while(fgets(Name_Buffer, bufLen, x->fp))
{
x->query_num++;
}
x->read_name = (char**)malloc(sizeof(char*)*x->query_num);
x->candidate_count = (kvec_t_u64_warp*)malloc(sizeof(kvec_t_u64_warp)*x->query_num);
fseek(x->fp, 0, SEEK_SET);
i = 0;
while(fgets(Name_Buffer, bufLen, x->fp))
{
nameLen = strlen(Name_Buffer) - 1;
x->read_name[i] = (char*)malloc(sizeof(char)*(nameLen+1));
memcpy(x->read_name[i], Name_Buffer, sizeof(char)*nameLen);
x->read_name[i][nameLen] = '\0';
kv_init(x->candidate_count[i].a);
i++;
}
fclose(x->fp);
sprintf(Name_Buffer, "%s.debug.stdout", file);
x->fp = fopen(Name_Buffer,"w");
fprintf(stderr, "Print debugging information to: %s\n", Name_Buffer);
free(Name_Buffer);
}
void destory_Debug_reads(Debug_reads* x)
{
uint64_t i;
for (i = 0; i < x->query_num; i++)
{
free(x->read_name[i]);
kv_destroy(x->candidate_count[i].a);
}
free(x->read_name);
fclose(x->fp);
}
+18 -2
View File
@@ -22,7 +22,8 @@
///#define Get_READ(R_INF, ID) R_INF.read + (R_INF.index[ID]>>2) + ID
#define Get_READ(R_INF, ID) (R_INF).read_sperate[(ID)]
#define Get_NAME(R_INF, ID) ((R_INF).name + (R_INF).name_index[(ID)])
#define CHECK_BY_NAME(R_INF, NAME, ID) (Get_NAME_LENGTH((R_INF),(ID))==strlen((NAME)) && \
memcmp((NAME), Get_NAME((R_INF), (ID)), Get_NAME_LENGTH((R_INF),(ID))) == 0)
extern uint8_t seq_nt6_table[256];
extern char bit_t_seq_table[256][4];
@@ -93,6 +94,7 @@ typedef struct
uint32_t new_length;
} Compressed_Cigar_record;
#define AMBIGU 0
#define FATHER 1
#define MOTHER 2
@@ -128,6 +130,8 @@ typedef struct
ma_hit_t_alloc* paf;
ma_hit_t_alloc* reverse_paf;
///kvec_t_u64_warp* pb_regions;
} All_reads;
extern All_reads R_INF;
@@ -140,6 +144,15 @@ typedef struct
long long RID;
} UC_Read;
typedef struct
{
char** read_name;
uint64_t query_num;
kvec_t_u64_warp* candidate_count;
FILE* fp;
pthread_mutex_t OutputMutex;
} Debug_reads;
void init_All_reads(All_reads* r);
void malloc_All_reads(All_reads* r);
void ha_insert_read_len(All_reads *r, int read_len, int name_len);
@@ -153,5 +166,8 @@ void reverse_complement(char* pattern, uint64_t length);
void write_All_reads(All_reads* r, char* read_file_name);
int load_All_reads(All_reads* r, char* read_file_name);
void destory_All_reads(All_reads* r);
int destory_read_bin(All_reads* r);
void init_Debug_reads(Debug_reads* x, const char* file);
void destory_Debug_reads(Debug_reads* x);
void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID);
#endif
+683 -17
View File
@@ -6,6 +6,10 @@
#include "Overlaps.h"
#include "Correct.h"
#include "kthread.h"
#include "kdq.h"
#include "hic.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)))
@@ -100,6 +104,7 @@ typedef struct {
ma_ug_t *ug;
asg_t *read_g;
ma_hit_t_alloc* sources;
ma_hit_t_alloc* reverse_sources;
R_to_U* ruIndex;
ma_sub_t *coverage_cut;
@@ -109,9 +114,352 @@ typedef struct {
int min_ovlp;
float chain_rate;
hap_overlaps_list* all_ovlp;
long long cov_threshold;
}hap_alignment_struct_pip;
void print_peak_line(int c, int x, int exceed, int64_t cnt)
{
int j;
if (c >= 0) fprintf(stderr, "[M::%s] %5d: ", __func__, c);
else fprintf(stderr, "[M::%s] %5s: ", __func__, "rest");
for (j = 0; j < x; ++j) fputc('*', stderr);
if (exceed) fputc('>', stderr);
fprintf(stderr, " %lld\n", (long long)cnt);
}
void print_peak(long long* cov_buf, long long cov_buf_length, long long max_i)
{
long long i;
const long long hist_max = 100;
// print histogram
for (i = 0; i < cov_buf_length; ++i)
{
long long x, exceed = 0;
x = (int)((double)hist_max * cov_buf[i] / cov_buf[max_i] + .499);
if (x > hist_max) exceed = 1, x = hist_max; // may happen if cnt[2] is higher
if (i > max_i && x == 0) break;
print_peak_line(i, x, exceed, cov_buf[i]);
}
{
long long x, exceed = 0;
long long rest = 0;
for (; i < cov_buf_length; ++i) rest += cov_buf[i];
x = (int)((double)hist_max * rest / cov_buf[max_i] + .499);
if (x > hist_max) exceed = 1, x = hist_max;
print_peak_line(-1, x, exceed, rest);
}
}
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;
i = start = err_i = max_i = max2_i = max3_i = topo_peak_i = -1;
max = max2 = max3 = topo_peak = min = -1;
///cov_buf[0] is usually very large
for (i = 1; i < cov_buf_length; ++i)
{
if(cov_buf[i] > cov_buf[i-1]) break;
}
err_i = i - 1;
// find the global highest peak
max_i = err_i + 1, max = cov_buf[max_i];
for (i = max_i; i < cov_buf_length; ++i)
{
if (cov_buf[i] > max)
{
max = cov_buf[i];
max_i = i;
}
}
///print_peak(cov_buf, cov_buf_length, max_i);
// look for smaller peak on the low end
max2 = -1; max2_i = -1;
for (i = max_i - 1; i > err_i; --i)
{
///at first, it should be a peak
if (cov_buf[i] >= cov_buf[i-1] && cov_buf[i] >= cov_buf[i+1])
{
if (cov_buf[i] > max2)
{
max2 = cov_buf[i];
max2_i = i;
}
}
}
///fprintf(stderr, "***max2: %lld, max2_i: %lld\n", max2, max2_i);
if (max2_i != -1 && max2_i > err_i && max2_i < max_i)
{
for (i = max2_i + 1, min = max; i < max_i; ++i)
{
if (cov_buf[i] < min) min = cov_buf[i];
}
///if the second peak is not significant
if(max2 < max * 0.05 || min > max2 * 0.95) max2 = max2_i = -1;
}
// look for smaller peak on the high end
max3 = -1; max3_i = -1;
// we'd better use i < cov_buf_length - 1, since cov_buf[cov_buf_length-1] may have problem
for (i = max_i + 1; i < cov_buf_length - 1; ++i)
{
//at first, it should be a peak
if (cov_buf[i] >= cov_buf[i-1] && cov_buf[i] >= cov_buf[i+1])
{
if (cov_buf[i] > max3)
{
max3 = cov_buf[i], max3_i = i;
}
}
}
///fprintf(stderr, "***max3: %lld, max3_i: %lld\n", max3, max3_i);
//if found a peak
if (max3 != -1 && max3_i > max_i)
{
for (i = max_i + 1, min = max; i < max3_i; ++i)
{
if (cov_buf[i] < min) min = cov_buf[i];
}
if (max3 < max * 0.05 || min > max3 * 0.95 || max3_i > max_i * 3) max3 = max3_i = -1;
}
(*hom_peak) = (*het_peak) = -1;
if (topo_peak_cov && (*topo_peak_cov) < cov_buf_length)
{
topo_peak_i = (*topo_peak_cov);
topo_peak = cov_buf[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;
alter_peak = topo_peak_i;
k_mer_het = asm_opt.het_cov;
k_mer_hom = asm_opt.hom_cov;
if(max3_i > 0)
{
coverage_het = max_i;
coverage_hom = max3_i;
}
else
{
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);
}
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)
{
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;
uint32_t tn, is_Unitig;
long long* cov_buf = NULL;
ma_hit_t *h;
cov_buf = (long long*)calloc(cov_buf_length, sizeof(long long));
long long R_bases = 0, C_bases_primary = 0, C_bases_alter = 0, C_bases = 0;
memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);
for (v = 0; v < n_vtx; ++v)
{
if(nsg->seq[v].del) continue;
if(nsg->seq[v].c == ALTER_LABLE) continue;
u = &(ug->u.a[v]);
if(u->m == 0) continue;
for (k = 0; k < u->n; k++)
{
qn = u->a[k]>>33;
position_index[qn] = 0;
R_bases = coverage_cut[qn].e - coverage_cut[qn].s;
primary_bases += R_bases;
}
}
for (qn = 0; qn < read_g->n_seq; qn++)
{
if(position_index[qn] == 0) continue;
if(read_g->seq[qn].del) continue;
C_bases = C_bases_primary = C_bases_alter = 0;
R_bases = coverage_cut[qn].e - coverage_cut[qn].s;
alter_bases += R_bases;
for (j = 0; j < (uint64_t)(sources[qn].length); j++)
{
h = &(sources[qn].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(position_index[tn] == 0)
{
C_bases_primary += Get_qe((*h)) - Get_qs((*h));
}
else
{
C_bases_alter += Get_qe((*h)) - Get_qs((*h));
}
}
// 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;
C_bases = C_bases/R_bases;
if(C_bases < 0 || C_bases >= cov_buf_length) continue;
cov_buf[C_bases]++;
}
long long max_i = -1, max = -1;
for (j = 0; (long long)j < cov_buf_length; ++j)
{
if (cov_buf[j] > max)
{
max = cov_buf[j];
max_i = j;
}
}
///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);
memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);
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* k_mer_only, long long* coverage_only)
{
uint64_t i, j;
long long* cov_buf = NULL;
ma_hit_t *h;
cov_buf = (long long*)calloc(cov_buf_length, sizeof(long long));
long long R_bases = 0, C_bases = 0;
for (i = 0; i < n_read; ++i)
{
C_bases = 0;
R_bases = coverage_cut[i].e - coverage_cut[i].s;
for (j = 0; j < (uint64_t)(sources[i].length); j++)
{
h = &(sources[i].buffer[j]);
if(h->el != 1) continue;
C_bases += Get_qe((*h)) - Get_qs((*h));
}
C_bases = C_bases/R_bases;
if(C_bases < 0 || C_bases >= cov_buf_length) continue;
cov_buf[C_bases]++;
}
long long alter_peak = -1, hom_peak = -1, het_peak = -1;
if(position_index)
{
alter_peak = get_alter_peak(ug, read_g, ruIndex, position_index, sources, coverage_cut,
cov_buf_length);
}
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*HOM_PEAK_RATE;
if(het_peak != -1) return het_peak*HET_PEAK_RATE;
return -1;
}
void init_hap_alignment_struct(hap_alignment_struct* x, uint32_t size)
{
x->vote_counting = (uint64_t*)malloc(sizeof(uint64_t)*size);
@@ -143,7 +491,7 @@ void destory_hap_alignment_struct(hap_alignment_struct* x)
}
void init_hap_alignment_struct_pip(hap_alignment_struct_pip* x, uint32_t num_threads, uint32_t n_seq,
ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, ma_sub_t *coverage_cut,
ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, ma_sub_t *coverage_cut,
uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chain_rate, hap_overlaps_list* all_ovlp)
{
uint32_t i;
@@ -156,6 +504,7 @@ uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chai
x->ug = ug;
x->read_g = read_g;
x->sources = sources;
x->reverse_sources = reverse_sources;
x->ruIndex = ruIndex;
x->coverage_cut = coverage_cut;
@@ -487,6 +836,27 @@ long long* r_yBeg, long long* r_yEnd)
return XCY;
}
uint64_t get_pair_hap_coverage(uint64_t* readIDs, uint32_t Len, ma_hit_t_alloc* sources, ma_sub_t* coverage_cut)
{
uint32_t m, n, qn;
ma_hit_t *h;
uint64_t R_bases = 0, C_bases = 0;
for (m = 0; m < Len; m++)
{
qn = readIDs[m]>>33;
R_bases += coverage_cut[qn].e - coverage_cut[qn].s;
for (n = 0; n < (uint64_t)(sources[qn].length); n++)
{
h = &(sources[qn].buffer[n]);
C_bases += Get_qe((*h)) - Get_qs((*h));
}
}
return C_bases/R_bases;
}
void get_pair_hap_similarity(uint64_t* readIDs, uint32_t Len, uint32_t target_uId,
ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match, double* Total)
{
@@ -1458,11 +1828,10 @@ long long* r_y_pos_beg, long long* r_y_pos_end)
uint32_t calculate_pair_hap_similarity_advance(hap_candidates* hap_can,
uint64_t* position_index, uint32_t xUid, uint32_t yUid, ma_utg_t* xReads, ma_utg_t* yReads,
ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, ma_sub_t *coverage_cut,
float Hap_rate, int max_hang, int min_ovlp, kvec_asg_arc_t_offset* u_buffer,
kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg,
long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end)
uint64_t* position_index, uint32_t xUid, uint32_t yUid, ma_utg_t* xReads, ma_utg_t* yReads,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, ma_sub_t *coverage_cut,
float Hap_rate, int max_hang, int min_ovlp, uint64_t cov_threshold, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex,
long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end)
{
uint32_t max_count = 0, min_count = 0, flag;
uint32_t xLen = xReads->n, xIndex;
@@ -1532,6 +1901,7 @@ long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end)
if(max_count > min_count*Hap_rate)
{
long long r_x_interval_beg, r_x_interval_end, r_y_interval_beg, r_y_interval_end;
uint64_t ploid_coverage = 0;
///for containment, don't need to do anything
get_hap_alignment_boundary(xReads, yReads, flag, xLeftMatch, xLeftTotal,
@@ -1565,6 +1935,15 @@ long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end)
if(hap_can->index_end == YCX && xReads->len > (yReads->len*2)) return NON_PLOID;
if(hap_can->index_end == (uint32_t)-1) return NON_PLOID;
ploid_coverage = 0;
ploid_coverage += get_pair_hap_coverage(xReads->a+r_x_interval_beg, r_x_interval_end+1-r_x_interval_beg,
sources, coverage_cut);
ploid_coverage += get_pair_hap_coverage(yReads->a+r_y_interval_beg, r_y_interval_end+1-r_y_interval_beg,
sources, coverage_cut);
///fprintf(stderr, "ploid_coverage: %lu, cov_threshold: %lu\n", ploid_coverage, cov_threshold);
if(cov_threshold > 0 && ploid_coverage >= cov_threshold) return NON_PLOID;
return PLOID;
}
return NON_PLOID;
@@ -2425,6 +2804,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid)
hap_alignment_struct_pip* hap_buf = (hap_alignment_struct_pip*)_data;
ma_ug_t *ug = hap_buf->ug;
asg_t *read_g = hap_buf->read_g;
ma_hit_t_alloc* sources = hap_buf->sources;
ma_hit_t_alloc* reverse_sources = hap_buf->reverse_sources;
R_to_U* ruIndex = hap_buf->ruIndex;
ma_sub_t *coverage_cut = hap_buf->coverage_cut;
@@ -2444,6 +2824,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);
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;
@@ -2591,9 +2973,9 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid)
if(u_can->a.a[k].weight < Get_match(hap_can)*Hap_rate) continue;
if(calculate_pair_hap_similarity_advance(&(u_can->a.a[k]), position_index, xUid, yUid,
xReads, yReads, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, max_hang,
min_ovlp, u_buffer, score_vc, prevIndex_vec, &r_x_pos_beg, &r_x_pos_end, &r_y_pos_beg,
&r_y_pos_end)!=PLOID)
xReads, yReads, sources, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, max_hang,
min_ovlp, cov_threshold, u_buffer, score_vc, prevIndex_vec, &r_x_pos_beg, &r_x_pos_end,
&r_y_pos_beg, &r_y_pos_end)!=PLOID)
{
continue;
}
@@ -2910,7 +3292,7 @@ int asg_pop_bubble_purge_graph(asg_t *purge_g, int max_dist)
for (i = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs
if (!av[i].del) ++n_arc;
if (n_arc > 1)
n_pop += asg_bub_pop1_primary_trio(purge_g, NULL, v, max_dist, &b, (uint32_t)-1, DROP, 1);
n_pop += asg_bub_pop1_primary_trio(purge_g, NULL, v, max_dist, &b, (uint32_t)-1, DROP, 1, NULL, NULL);
}
free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
if (n_pop) asg_cleanup(purge_g);
@@ -3544,10 +3926,57 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg
}
void collect_reverse_unitig_pair(hc_links* link, ma_ug_t *ug, hap_overlaps* t)
{
uint32_t i = 0, k = 0, rId_0, rId_1, pre_0, pre_1, b_0 = t->xUid, b_1 = t->yUid;
uint64_t d = RC_2;
ma_utg_t* u_b_0 = &(ug->u.a[b_0]);
ma_utg_t* u_b_1 = &(ug->u.a[b_1]);
if(u_b_0->n == 0) return;
if(u_b_1->n == 0) return;
for (i = t->x_beg_id, pre_0 = (uint32_t)-1; i < t->x_end_id; i++)
{
rId_0 = u_b_0->a[i]>>33;
if(link->u_idx[rId_0] == (uint32_t)-1) continue;
if(pre_0 == link->u_idx[rId_0]) continue;
pre_0 = link->u_idx[rId_0];
for (k = t->y_beg_id, pre_1 = (uint32_t)-1; k < t->y_end_id; k++)
{
rId_1 = u_b_1->a[k]>>33;
if(link->u_idx[rId_1] == (uint32_t)-1) continue;
if(pre_1 == link->u_idx[rId_1]) continue;
pre_1 = link->u_idx[rId_1];
push_hc_edge(&(link->a.a[pre_0]), pre_1, 1, 1, &d);
push_hc_edge(&(link->a.a[pre_1]), pre_0, 1, 1, &d);
}
}
}
void collect_reverse_unitigs_purge(buf_t* b_0, hc_links* link, ma_ug_t *ug, hap_overlaps_list* all_ovlp)
{
if(b_0->b.n <= 1) return;
uint32_t k;
int index = 0;
for (k = 0; k < b_0->b.n - 1; k++)
{
index = get_specific_hap_overlap(&(all_ovlp->x[b_0->b.a[k]>>1]), b_0->b.a[k]>>1, b_0->b.a[k+1]>>1);
if(index == -1) continue;
collect_reverse_unitig_pair(link, ug, &(all_ovlp->x[b_0->b.a[k]>>1].a.a[index]));
}
}
void link_unitigs(asg_t *purge_g, ma_ug_t *ug, hap_overlaps_list* all_ovlp,
R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g,
uint64_t* position_index, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex,
kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit)
kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit,
hc_links* link)
{
uint32_t v, n_vtx = purge_g->n_seq * 2, beg, end;
long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen;
@@ -3567,6 +3996,8 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg
{
continue;
}
if(link) collect_reverse_unitigs_purge(&b_0, link, ug, all_ovlp);
purge_merge(purge_g, ug, all_ovlp, &b_0, ruIndex, reverse_sources, coverage_cut,
read_g, position_index, u_buffer, tailIndex, prevIndex,max_hang, min_ovlp, edge, visit);
}
@@ -3587,9 +4018,210 @@ void print_all_purge_ovlp(ma_ug_t *ug, hap_overlaps_list* all_ovlp)
}
}
void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, 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, uint32_t just_contain)
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,
uint32_t just_contain, uint32_t just_coverage, hc_links* link)
{
asg_t *purge_g = NULL;
purge_g = asg_init();
@@ -3620,7 +4252,19 @@ int min_ovlp, long long bubble_dist, float drop_ratio, uint32_t just_contain)
asg_arc_t* p = NULL;
int r;
hap_alignment_struct_pip hap_buf;
long long k_mer_only, coverage_only;
if(asm_opt.hom_global_coverage != -1)
{
hap_buf.cov_threshold = asm_opt.hom_global_coverage;
}
else
{
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);
}
for (v = 0; v < nsg->n_seq; v++)
{
uId = v;
@@ -3634,7 +4278,7 @@ int min_ovlp, long long bubble_dist, float drop_ratio, uint32_t just_contain)
for (i = 0, offset = 0; i < reads->n; i++)
{
rId = reads->a[i]>>33;
set_R_to_U(ruIndex, rId, uId, 1);
set_R_to_U(ruIndex, rId, uId, 1, &(read_g->seq[rId].c));
position_index[rId] = offset;
position_index[rId] = position_index[rId] << 32;
@@ -3649,9 +4293,27 @@ int min_ovlp, long long bubble_dist, float drop_ratio, uint32_t just_contain)
init_hap_alignment_struct_pip(&hap_buf, asm_opt.thread_num, nsg->n_seq, ug, read_g,
reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp,
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;
}
}
if(asm_opt.hom_global_coverage == -1) asm_opt.hom_global_coverage = hap_buf.cov_threshold;
fprintf(stderr, "[M::%s] purge duplication coverage threshold: %lld\n", __func__, hap_buf.cov_threshold);
if(just_coverage) goto end_coverage;
///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);
@@ -3688,6 +4350,7 @@ int min_ovlp, long long bubble_dist, float drop_ratio, uint32_t just_contain)
purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].c = ALTER_LABLE;
purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].del = 1;
all_ovlp.x[uId].a.a[i].status = DELETE;
if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i]));
}
if(all_ovlp.x[uId].a.a[i].type == XCY)
@@ -3696,6 +4359,7 @@ int min_ovlp, long long bubble_dist, float drop_ratio, uint32_t just_contain)
purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].c = ALTER_LABLE;
purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].del = 1;
all_ovlp.x[uId].a.a[i].status = DELETE;
if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i]));
}
///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
}
@@ -3747,7 +4411,7 @@ int min_ovlp, long long bubble_dist, float drop_ratio, uint32_t just_contain)
link_unitigs(purge_g, ug, &all_ovlp, ruIndex, reverse_sources, coverage_cut, read_g, position_index,
&(hap_buf.buf[0].u_buffer), &(hap_buf.buf[0].u_buffer_tailIndex), &(hap_buf.buf[0].u_buffer_prevIndex),
max_hang, min_ovlp, edge, hap_buf.buf[0].visit);
max_hang, min_ovlp, edge, hap_buf.buf[0].visit, link);
}
for (v = 0; v < all_ovlp.num; v++)
@@ -3759,6 +4423,8 @@ int min_ovlp, long long bubble_dist, float drop_ratio, uint32_t just_contain)
}
}
end_coverage:
uint32_t is_Unitig;
for (v = 0; v < ruIndex->len; v++)
{
+10 -4
View File
@@ -6,13 +6,19 @@
#include "kdq.h"
#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* 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, uint32_t just_contain);
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,
uint32_t just_contain, uint32_t just_coverage, hc_links* link);
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
+144 -63
View File
@@ -4,18 +4,102 @@
# Install hifiasm (requiring g++ and zlib)
git clone https://github.com/chhylp123/hifiasm
cd hifiasm && make
# Assembly
./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz
# Run on test data (use -f0 for small datasets)
wget https://github.com/chhylp123/hifiasm/releases/download/v0.7/chr11-2M.fa.gz
./hifiasm -o test -t4 -f0 chr11-2M.fa.gz 2> test.log
awk '/^S/{print ">"$2;print $3}' test.p_ctg.gfa > test.p_ctg.fa # get primary contigs in FASTA
# Assemble inbred/homozygous genomes (-l0 disables duplication purging)
hifiasm -o CHM13.asm -t32 -l0 CHM13-HiFi.fa.gz 2> CHM13.asm.log
# Assemble heterozygous with built-in duplication purging
hifiasm -o HG002.asm -t32 HG002-file1.fq.gz HG002-file2.fq.gz
# Trio binning assembly (requiring https://github.com/lh3/yak)
yak count -b37 -t16 -o pat.yak <(cat pat_1.fq.gz pat_2.fq.gz) <(cat pat_1.fq.gz pat_2.fq.gz)
yak count -b37 -t16 -o mat.yak <(cat mat_1.fq.gz mat_2.fq.gz) <(cat mat_1.fq.gz mat_2.fq.gz)
hifiasm -o HG002.asm -t32 -1 pat.yak -2 mat.yak HG002-HiFi.fa.gz
```
## Introduction
Hifiasm is a fast haplotype-resolved de novo assembler for PacBio
Hifi reads. Unlike most existing assemblers, hifiasm starts from uncollapsed
genome. Thus, it is able to keep the haplotype information as much as possible.
Hifiasm is a fast haplotype-resolved de novo assembler for PacBio Hifi reads.
It can assemble a human genome in several hours and works with the California
redwood genome, one of the most complex genomes sequenced so far. Hifiasm can
produce primary/alternate assemblies of quality competitive with the best
assemblers. It also introduces a new graph binning algorithm and achieves
the best haplotype-resolved assembly given trio data.
For non-trio assembly, the input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and its
outputs consist of:
## Why Hifiasm?
* Hifiasm delivers high-quality assemblies. It tends to generate longer contigs
and resolve more segmental duplications than other assemblers.
* Given sequence reads from the parents, hifiasm can produce overall the best
haplotype-resolved assembly so far. It is the assembler of choice by the
[Human Pangenome Project][hpp] for the first batch of samples.
* Hifiasm can purge duplications between haplotigs without relying on
third-party tools such as purge\_dups. Hifiasm does not need polishing tools
like pilon or racon, either. This simplifies the assembly pipeline and saves
running time.
* Hifiasm is fast. It can assemble a human genome in half a day and assemble a
~30Gb redwood genome in three days. No genome is too large for hifiasm.
* Hifiasm is trivial to install and easy to use. It does not required python,
R or C++11 compilers and can be compiled into a single executable. The
default setting works well with a variety of genomes.
[hpp]: https://humanpangenome.org
## Usage
A typical hifiasm command line looks like:
```sh
hifiasm -o NA12878.asm -t 32 NA12878.fq.gz
```
where `NA12878.fq.gz` provides the input reads, `-t` sets the number of CPUs in
use and `-o` specifies the prefix of output files. For this example, the
primary contigs are written to `NA12878.asm.p_ctg.gfa` and alternate contigs to
`NA12878.asm.a_ctg.gfa`. At the first run, hifiasm saves corrected reads and
overlaps to disk as `NA12878.asm.*.bin`. It reuses the saved results to avoid
the time-consuming all-vs-all overlap calculation next time. You may specify
`-i` to ignore precomputed overlaps and redo overlapping from raw reads.
Hifiasm purges haplotig duplications by default. For inbred or homozygous
genomes, you may disable purging with option `-l0`. Old HiFi reads may contain
short adapter sequences at the ends of reads. You can specify `-z20` to trim
both ends of reads by 20bp. For small genomes, use `-f0` to disable the initial
bloom filter which takes 16GB memory at the beginning. For genomes much larger
than human, applying `-f38` or even `-f39` is preferred to save memory on k-mer
counting.
When parental short reads are available, hifiasm can generate a pair of
haplotype-resolved assemblies with trio binning. To perform such assembly, you
need to count k-mers first with [yak][yak] first and then do assembly:
```sh
yak count -k31 -b37 -t16 -o pat.yak paternal.fq.gz
yak count -k31 -b37 -t16 -o mat.yak maternal.fq.gz
hifiasm -o NA12878.asm -t 32 -1 pat.yak -2 mat.yak NA12878.fq.gz
```
Here `NA12878.asm.hap1.p_ctg.gfa` and `NA12878.asm.hap2.p_ctg.gfa` give the two
haplotype assemblies. In the binning mode, hifiasm does not purge haplotig
duplications by default. Because hifiasm reuses saved overlaps, you can
generate both primary/alternate assemblies and trio binning assemblies with
```sh
hifiasm -o NA12878.asm -t 32 NA12878.fq.gz 2> NA12878.asm.pri.log
hifiasm -o NA12878.asm -t 32 -1 pat.yak -2 mat.yak /dev/null 2> NA12878.asm.trio.log
```
The second command line will run much faster than the first. You can also dump
error corrected in FASTA and/or overlaps in PAF with
```sh
hifiasm -o NA12878.asm -t 32 --write-paf --write-ec /dev/null
```
## Output files
For non-trio assembly, hifiasm generates the following files:
1. Haplotype-resolved raw [unitig][unitig] graph in [GFA][gfa] format
(*prefix*.r\_utg.gfa). This graph keeps all haplotype information, including
@@ -28,7 +112,8 @@ outputs consist of:
4. Alternate assembly contig graph (*prefix*.a\_ctg.gfa). This graph consists of all assemblies that
are discarded in primary contig graph.
For trio assembly, the input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and the paternal/maternal trio indexes generated by `yak count` (see https://github.com/lh3/yak). The outputs consist of:
For trio assembly, hifiasm generates the following files:
1. Haplotype-resolved raw [unitig][unitig] graph in [GFA][gfa] format
(*prefix*.r\_utg.gfa). This graph keeps all haplotype information.
@@ -38,77 +123,73 @@ For trio assembly, the input of hifiasm is the PacBio Hifi reads in fasta/fastq
3. Phased maternal/haplotype2 contig graph (*prefix*.hap2.p\_ctg.gfa). This graph keeps the phased
maternal/haplotype2 assembly.
Hifiasm writes error corrected reads to the *prefix*.ec.bin binary file and
writes overlaps to *prefix*.ovlp.source.bin and *prefix*.ovlp.reverse.bin.
## Results
In addition, hifiasm also outputs three binary files that save all overlap information (*prefix*.ec.bin, *prefix*.ovlp.reverse.bin, *prefix*.ovlp.source.bin). With these files, hifiasm can avoid the time-consuming all-to-all overlap calculation step, and do the assembly
directly and quickly. This might be helpful when you want to get an optimized
assembly by multiple rounds of experiments with different parameters.
The following table shows the statistics of several hifiasm primary assemblies:
Hifiasm is a standalone and lightweight assembler, which does not need external
libraries (except zlib). For large genomes, it can generate high-quality
assembly in a few hours. Hifiasm has been tested on human, butterfly, rice and drosophila.
In particular, hifiasm is able to assemble the 26.5Gb California redwood tree in a few days.
The results are as follows:
|<sub>Dataset<sub>|<sub>GSize<sub>|<sub>Cov<sub>|<sub>Asm options<sub>|<sub>CPU time<sub>|<sub>Wall time<sub>|<sub>RAM<sub>|<sub>[unitig][unitig]/[contig][unitig] N50<sup>[1]</sup><sub>|
|<sub>Dataset<sub>|<sub>Size<sub>|<sub>Cov.<sub>|<sub>Asm options<sub>|<sub>CPU time<sub>|<sub>Wall time<sub>|<sub>RAM<sub>|<sub> N50<sub>|
|:---------------|-----:|-----:|:---------------------|-------:|--------:|----:|----------------:|
|<sub>[\[Redwood\]](https://downloads.pacbcloud.com/public/dataset/redwood2020/)<sub>|<sub>26.5Gb<sub>|<sub>x23<sub>|<sub>-k 40 -t 64 -r 2<sub>|<sub>7274h30m<sub>|<sub>141h30m<sub>|<sub>512G<sub>|<sub>1.7Mb/1.9Mb<sub>|
|<sub>[Mouse (C57/BL6J)][mouse-data]</sub>|<sub>2.6Gb</sub> |<sub>&times;25</sub>|<sub>-t48 -l0</sub> |<sub>172.9h</sub> |<sub>4.8h</sub> |<sub>76G</sub> |<sub>21.1Mb</sub>|
|<sub>[Maize (B73)][maize-data]</sub> |<sub>2.2Gb</sub> |<sub>&times;22</sub>|<sub>-t48 -l0</sub> |<sub>203.2h</sub> |<sub>5.1h</sub> |<sub>68G</sub> |<sub>36.7Mb</sub>|
|<sub>[Strawberry][strawberry-data]</sub> |<sub>0.8Gb</sub> |<sub>&times;36</sub>|<sub>-t48 -D10</sub>|<sub>152.7h</sub> |<sub>3.7h</sub> |<sub>91G</sub> |<sub>17.8Mb</sub>|
|<sub>[Frog][frog-data]</sub> |<sub>9.5Gb</sub> |<sub>&times;29</sub>|<sub>-t48</sub> |<sub>2834.3h</sub>|<sub>69.0h</sub>|<sub>463G</sub>|<sub>9.3Mb</sub>|
|<sub>[Redwood][redwood-data]</sub> |<sub>35.6Gb</sub>|<sub>&times;28</sub>|<sub>-t80</sub> |<sub>3890.3h</sub>|<sub>65.5h</sub>|<sub>699G</sub>|<sub>5.4Mb</sub>|
|<sub>[Human (CHM13)][CHM13-data]</sub> |<sub>3.1Gb</sub> |<sub>&times;32</sub>|<sub>-t48 -l0</sub> |<sub>310.7h</sub> |<sub>8.2h</sub> |<sub>114G</sub>|<sub>88.9Mb</sub>|
|<sub>[Human (HG00733)][HG00733-data]</sub>|<sub>3.1Gb</sub>|<sub>&times;33</sub>|<sub>-t48</sub> |<sub>269.1h</sub> |<sub>6.9h</sub> |<sub>135G</sub>|<sub>69.9Mb</sub>|
|<sub>[Human (HG002)][NA24385-data]</sub> |<sub>3.1Gb</sub> |<sub>&times;36</sub>|<sub>-t48</sub> |<sub>305.4h</sub> |<sub>7.7h</sub> |<sub>137G</sub>|<sub>98.7Mb</sub>|
<sub>[1] unitig N50 is the N50 of assembly graph with haplotype information (i.e., bubbles), while the contig N50 is the N50 of haplotype collapsed assembly (i.e., without bubbles).<sub>
[mouse-data]: https://www.ncbi.nlm.nih.gov/sra/?term=SRR11606870
[maize-data]: https://www.ncbi.nlm.nih.gov/sra/?term=SRR11606869
[strawberry-data]: https://www.ncbi.nlm.nih.gov/sra/?term=SRR11606867
[frog-data]: https://www.ncbi.nlm.nih.gov/sra?term=(SRR11606868)%20OR%20SRR12048570
[redwood-data]: https://www.ncbi.nlm.nih.gov/sra/?term=SRP251156
[CHM13-data]: https://www.ncbi.nlm.nih.gov/sra?term=(((SRR11292120)%20OR%20SRR11292121)%20OR%20SRR11292122)%20OR%20SRR11292123
## Usage
Hifiasm can assemble a 3.1Gb human genome in several hours or a ~30Gb hexaploid
redwood genome in a few days on a single machine. For trio binning assembly:
For Hifi reads assembly, a typical command line looks like:
|<sub>Dataset<sub>|<sub>Cov.<sub>|<sub>CPU time<sub>|<sub>Elapsed time<sub>|<sub>RAM<sub>|<sub> N50<sub>|
|:---------------|-----:|-------:|--------:|----:|----------------:|
|<sub>[HG00733][HG00733-data], [\[father\]][HG00731-data], [\[mother\]][HG00732-data]</sub>|<sub>&times;33</sub>|<sub>269.1h</sub>|<sub>6.9h</sub>|<sub>135G</sub>|<sub>35.1Mb (paternal), 34.9Mb (maternal)</sub>|
|<sub>[HG002][NA24385-data], [\[father\]][NA24149-data], [\[mother\]][NA24143-data]</sup>|<sub>&times;36</sub>|<sub>305.4h</sub>|<sub>7.7h</sub>|<sub>137G</sub>|<sub>41.0Mb (paternal), 40.8Mb (maternal)</sub>|
|<sub>[NA12878][NA12878-data], [\[father\]][NA12891-data], [\[mother\]][NA12892-data]</sub>|<sub>&times;30</sub>|<sub>180.8h</sub>|<sub>4.9h</sub>|<sub>123G</sub>|<sub>27.7Mb (paternal), 27.0Mb (maternal)</sub>|
```sh
./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz
[HG00733-data]: https://www.ebi.ac.uk/ena/data/view/ERX3831682
[HG00731-data]: https://www.ebi.ac.uk/ena/data/view/ERR3241754
[HG00732-data]: https://www.ebi.ac.uk/ena/data/view/ERR3241755
[NA24385-data]: https://www.ncbi.nlm.nih.gov/sra?term=(((SRR10382244)%20OR%20SRR10382245)%20OR%20SRR10382248)%20OR%20SRR10382249
[NA24149-data]: https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/HG003_NA24149_father/NIST_HiSeq_HG003_Homogeneity-12389378/HG003Run01-13262252/
[NA24143-data]: https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/HG004_NA24143_mother/NIST_HiSeq_HG004_Homogeneity-14572558/HG004Run01-15133132/
[NA12878-data]: https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/NA12878/PacBio_SequelII_CCS_11kb/
[NA12891-data]: https://www.ebi.ac.uk/ena/data/view/ERR194160
[NA12892-data]: https://www.ebi.ac.uk/ena/data/view/ERR194161
Except NA12878, the assemblies above were produced by hifiasm v0.12 and can be
downloaded at
```txt
ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/submission/hifiasm-0.12/
```
NA12878 was assembled with an older version of hifiasm and is available at
```txt
ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/NA12878-r253/
```
where `NA12878.fq.gz` is the input reads and `-o` specifies the output files.
In this example, all output files can be found at `NA12878.asm.*`. `-t` specifies
the number of CPU threads. Note that at first run, hifiasm will save all overlaps
to disk, which can avoid the time-consuming all-to-all overlap calculation next time.
For hifiasm, once the overlap information has been obtained during the previous run
in advance, it is able to load all overlaps from disk and then directly do assembly.
If you want to ignore the pre-computed overlap information, please specify `-i`.
Please note that some old Hifi reads may consist of short adapters. To improve
the assembly quality, adapters should be removed by `-z` as follow:
```sh
./hifiasm -o butterfly.asm -t 42 -z 20 butterfly.fq.gz
```
In this example, hifiasm will remove 20 bases from both ends of each read.
For trio assembly, first the trio indexes of paternal/maternal should be generated by
`yak count` (see https://github.com/lh3/yak):
```sh
./yak count -k31 -b37 -t16 -o mat.yak mat.fq.gz
```
```sh
./yak count -k31 -b37 -t16 -o pat.yak pat.fq.gz
```
and then run hifiasm as follow:
```sh
./hifiasm -o NA12878.asm -t 32 -1 pat.yak -2 mat.yak NA12878_1.fq.gz NA12878_2.fq.gz
```
[unitig]: http://wgs-assembler.sourceforge.net/wiki/index.php/Celera_Assembler_Terminology
[gfa]: https://github.com/pmelsted/GFA-spec/blob/master/GFA-spec.md
[paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
[yak]: https://github.com/lh3/yak
## Getting Help
For detailed description of options, please see `man ./hifiasm.1`.
The `-h` option of hifiasm also provides simple description of options. If you
have further questions, please raise an issue at the issue page.
For detailed description of options, please see `man ./hifiasm.1`. The `-h`
option of hifiasm also provides brief description of options. If you have
further questions, please raise an issue at the [issue
page](https://github.com/chhylp123/hifiasm/issues).
## Limitations and future works
## Limitations
1. The running time and memory usage should be further reduced.
2. The N50 should be further improved.
1. Purging haplotig duplications may introduce misassemblies.
+8 -5
View File
@@ -67,7 +67,7 @@ static yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ..
{
va_list ap;
FILE *fp;
uint32_t t[3];
uint32_t t[3], f_tmp = 0;
char magic[4];
int i, j, absent, min_cnt = 0, mid_cnt = 0, mode_err = 0;
uint64_t mask = (1ULL<<YAK_COUNTER_BITS) - 1, n_ins = 0, n_new = 0;
@@ -92,33 +92,36 @@ static yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ..
fclose(fp);
return 0;
}
fread(t, 4, 3, fp);
f_tmp += fread(t, 4, 3, fp);
if (t[2] != YAK_COUNTER_BITS) {
fprintf(stderr, "ERROR: saved counter bits: %d; compile-time counter bits: %d\n", t[2], YAK_COUNTER_BITS);
fclose(fp);
return 0;
}
///t[0] = k; t[1] = pre, t[2] = YAK_COUNTER_BITS;
ch = ch0 == 0? yak_ch_init(t[0], t[1]) : ch0;
assert((int)t[0] == ch->k && (int)t[1] == ch->pre);
for (i = 0; i < 1<<ch->pre; ++i) {
yak_ht_t *h = ch->h[i].h;
fread(t, 4, 2, fp);
f_tmp += fread(t, 4, 2, fp);
///t[0] = kh_capacity(h), t[1] = kh_size(h);
if (ch0 == 0) yak_ht_resize(h, t[0]);
for (j = 0; j < (int)t[1]; ++j) {
uint64_t key;
fread(&key, 8, 1, fp);
f_tmp += fread(&key, 8, 1, fp);
if (mode == YAK_LOAD_ALL) {
++n_ins;
yak_ht_put(h, key, &absent);
if (absent) ++n_new;
} else if (mode == YAK_LOAD_TRIOBIN1 || mode == YAK_LOAD_TRIOBIN2) {
int cnt = key & mask, x, shift = mode == YAK_LOAD_TRIOBIN1? 0 : 2;
//1. filter singleton k-mer; 2. label non-repeat and repeat
if (cnt >= mid_cnt) x = 2<<shift;
else if (cnt >= min_cnt) x = 1<<shift;
else x = -1;
if (x >= 0) {
khint_t k;
///no need cnt at all
key = (key & ~mask) | x;
++n_ins;
k = yak_ht_put(h, key, &absent);
+158 -20
View File
@@ -28,8 +28,8 @@ typedef struct {
} seed1_t;
struct ha_abuf_s {
uint64_t n_a, m_a;
uint32_t old_mz_m;
uint64_t n_a, m_a;///number of anchors (seed positions)
uint32_t old_mz_m;///number of seeds
ha_mz1_v mz;
seed1_t *seed;
anchor1_t *a;
@@ -50,11 +50,17 @@ uint64_t ha_abuf_mem(const ha_abuf_t *ab)
return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(ha_mz1_t) + sizeof(seed1_t)) + sizeof(ha_abuf_t);
}
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)
int ha_ov_type(const overlap_region *r, uint32_t len)
{
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
uint32_t i;
if (r->x_pos_s == 0 && r->x_pos_e == len - 1) return 2; // contained in a longer read
else if (r->x_pos_s > 0 && r->x_pos_e < len - 1) return 3; // containing a shorter read
else return r->x_pos_s == 0? 0 : 1;
}
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)
{
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);
@@ -64,9 +70,11 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
clear_overlap_region_alloc(overlap_list);
recover_UC_Read(ucr, &R_INF, rid);
ab->mz.n = 0, ab->n_a = 0;
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);
ha_sketch_query(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, k_flag, dbg_ct);
// minimizer of queried read
if (ab->mz.m > ab->old_mz_m) {
ab->old_mz_m = ab->mz.m;
REALLOC(ab->seed, ab->old_mz_m);
@@ -85,6 +93,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
}
for (i = 0, k = 0; i < ab->mz.n; ++i) {
int j;
///z is one of the minimizer
ha_mz1_t *z = &ab->mz.a[i];
seed1_t *s = &ab->seed[i];
for (j = 0; j < s->n; ++j) {
@@ -108,6 +117,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
}
}
// copy over to _cl_
if (ab->m_a >= (uint64_t)cl->size) {
cl->size = ab->m_a;
@@ -123,33 +133,34 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
}
cl->length = ab->n_a;
calculate_overlap_region_by_chaining(cl, overlap_list, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain);
calculate_overlap_region_by_chaining(cl, overlap_list, chain_idx, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain, f_cigar);
#if 0
if (overlap_list->length > 2000) {
fprintf(stderr, "B\t%ld\t%ld\t%ld\n", (long)rid, (long)overlap_list->length, (long)Get_READ_LENGTH(R_INF, rid));
if (overlap_list->length > 0) {
fprintf(stderr, "B\t%ld\t%ld\t%d\n", (long)rid, (long)overlap_list->length, rlen);
for (int i = 0; i < (int)overlap_list->length; ++i) {
overlap_region *r = &overlap_list->list[i];
fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand],
(int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed);
fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand],
(int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, ha_ov_type(r, rlen));
}
}
#endif
if ((int)overlap_list->length > max_n_chain) {
int32_t n[2], s[2];
n[0] = n[1] = 0, s[0] = s[1] = 0;
int32_t w, n[4], s[4];
n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0;
ks_introsort_or_ss(overlap_list->length, overlap_list->list);
for (i = 0; i < (uint32_t)overlap_list->length; ++i) {
const overlap_region *r = &overlap_list->list[i];
int dir = r->x_pos_s == 0? 0 : 1;
++n[dir];
if ((int)n[dir] == max_n_chain) s[dir] = r->shared_seed;
w = ha_ov_type(r, rlen);
++n[w];
if ((int)n[w] == max_n_chain) s[w] = r->shared_seed;
}
if (s[0] > 0 || s[1] > 0) {
if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) {
overlap_region *r = &overlap_list->list[i];
int dir = r->x_pos_s == 0? 0 : 1;
if (r->shared_seed > s[dir]) {
w = ha_ov_type(r, rlen);
if (r->shared_seed >= s[w]) {
if ((uint32_t)k != i) {
overlap_region t;
t = overlap_list->list[k];
@@ -163,5 +174,132 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
}
}
///ks_introsort_or_xs(overlap_list->length, overlap_list->list);
}
void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf)
{
uint64_t j = 0, inner_j = 0;
while (j < overlap_list->length && inner_j < paf->length)
{
if(overlap_list->list[j].y_id < paf->buffer[inner_j].tn)
{
j++;
}
else if(overlap_list->list[j].y_id > paf->buffer[inner_j].tn)
{
inner_j++;
}
else
{
if(overlap_list->list[j].y_pos_strand == paf->buffer[inner_j].rev)
{
overlap_list->list[j].is_match = 1;
}
j++;
inner_j++;
}
}
}
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)
{
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);
if(ha_idx_hp)
{
uint32_t i, k, y_id, overlapLen, max_i;
int shared_seed;
overlap_region t;
overlap_region_sort_y_id(overlap_list->list, overlap_list->length);
ma_hit_sort_tn(paf->buffer, paf->length);
ma_hit_sort_tn(rev_paf->buffer, rev_paf->length);
lable_matched_ovlp(overlap_list, paf);
lable_matched_ovlp(overlap_list, rev_paf);
for (i = 0, k = 0; i < overlap_list->length; ++i)
{
if(overlap_list->list[i].is_match == 1)
{
if(k != i)
{
t = overlap_list->list[k];
overlap_list->list[k] = overlap_list->list[i];
overlap_list->list[i] = t;
overlap_list->list[k].is_match = 0;
}
k++;
}
}
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);
if(overlap_list->length + overlap_list_hp->length > overlap_list->size)
{
overlap_list->list = (overlap_region*)realloc(overlap_list->list,
sizeof(overlap_region)*(overlap_list->length + overlap_list_hp->length));
memset(overlap_list->list + overlap_list->size, 0, sizeof(overlap_region)*
(overlap_list->length + overlap_list_hp->length - overlap_list->size));
overlap_list->size = overlap_list->length + overlap_list_hp->length;
}
for (i = 0, k = overlap_list->length; i < overlap_list_hp->length; i++, k++)
{
t = overlap_list->list[k];
overlap_list->list[k] = overlap_list_hp->list[i];
overlap_list_hp->list[i] = t;
}
overlap_list->length = k;
overlap_region_sort_y_id(overlap_list->list, overlap_list->length);
i = k = 0;
while (i < overlap_list->length)
{
y_id = overlap_list->list[i].y_id;
shared_seed = overlap_list->list[i].shared_seed;
overlapLen = overlap_list->list[i].overlapLen;
max_i = i;
i++;
while (i < overlap_list->length && overlap_list->list[i].y_id == y_id)
{
if((overlap_list->list[i].shared_seed > shared_seed) ||
((overlap_list->list[i].shared_seed == shared_seed) && (overlap_list->list[i].overlapLen <= overlapLen)))
{
y_id = overlap_list->list[i].y_id;
shared_seed = overlap_list->list[i].shared_seed;
overlapLen = overlap_list->list[i].overlapLen;
max_i = i;
}
i++;
}
if(k != max_i)
{
t = overlap_list->list[k];
overlap_list->list[k] = overlap_list->list[max_i];
overlap_list->list[max_i] = t;
}
k++;
}
overlap_list->length = k;
}
ks_introsort_or_xs(overlap_list->length, overlap_list->list);
}
void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list)
{
ks_introsort_or_xs(overlap_list->length, overlap_list->list);
}
+173
View File
@@ -0,0 +1,173 @@
#include <zlib.h>
#include <string.h>
#include "Process_Read.h"
#include "khashl.h"
#include "kseq.h"
typedef const char *cstr_t;
KHASHL_CSET_INIT(KH_LOCAL, strset_t, ss, cstr_t, kh_hash_str, kh_eq_str)
KHASHL_MAP_INIT(KH_LOCAL, hm64_t, h64, uint64_t, int, kh_hash_uint64, kh_eq_generic)
KSTREAM_INIT(gzFile, gzread, 65536)
#define GFA_MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr))))
#define GFA_REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr))))
char *gfa_strdup(const char *src)
{
int32_t len;
char *dst;
len = strlen(src);
GFA_MALLOC(dst, len + 1);
memcpy(dst, src, len + 1);
return dst;
}
char *gfa_strndup(const char *src, size_t n)
{
char *dst;
GFA_MALLOC(dst, n + 1);
strncpy(dst, src, n);
dst[n] = 0;
return dst;
}
char **gv_read_list(const char *o, int *n_)
{
int n = 0, m = 0;
char **s = 0;
*n_ = 0;
if (*o != '@') {
const char *q = o, *p;
for (p = q;; ++p) {
if (*p == ',' || *p == 0) {
if (n == m) {
m = m? m<<1 : 16;
GFA_REALLOC(s, m);
}
s[n++] = gfa_strndup(q, p - q);
if (*p == 0) break;
q = p + 1;
}
}
} else {
gzFile fp;
kstream_t *ks;
kstring_t str = {0,0,0};
int dret;
fp = gzopen(o + 1, "r");
if (fp == 0) return 0;
ks = ks_init(fp);
while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) {
char *p;
for (p = str.s; *p && !isspace(*p); ++p);
if (n == m) {
m = m? m<<1 : 16;
GFA_REALLOC(s, m);
}
s[n++] = gfa_strndup(str.s, p - str.s);
}
ks_destroy(ks);
gzclose(fp);
}
if (s) s = (char**)realloc(s, n * sizeof(char*));
*n_ = n;
return s;
}
void ha_extract_print(const All_reads *rs, int n_rounds, int n, char **list)
{
hm64_t *h;
khint_t k;
int i, absent, m, l;
uint64_t j;
const ma_hit_t_alloc *ov[2] = { rs->paf, rs->reverse_paf };
FILE *fp = stdout;
if (n > 0) {
int max_len = 0;
char *s = 0;
strset_t *ss;
ss = ss_init();
for (i = 0; i < n; ++i)
ss_put(ss, list[i], &absent);
for (j = 0; j < rs->total_reads; ++j)
if (max_len < (int)Get_NAME_LENGTH(*rs, j))
max_len = Get_NAME_LENGTH(*rs, j);
GFA_MALLOC(s, max_len + 1);
h = h64_init();
for (j = 0; j < rs->total_reads; ++j) {
strncpy(s, Get_NAME(*rs, j), Get_NAME_LENGTH(*rs, j));
s[Get_NAME_LENGTH(*rs, j)] = 0;
if (ss_get(ss, s) != kh_end(ss)) {
k = h64_put(h, j, &absent);
kh_val(h, k) = -1;
}
}
free(s);
ss_destroy(ss);
} else return;
for (m = 0; m < n_rounds; ++m) {
for (j = 0; j < rs->total_reads; ++j) {
for (l = 0; l < 2; ++l) {
const ma_hit_t_alloc *o = &ov[l][j];
for (i = 0; i < (int)o->length; ++i) {
uint64_t q = Get_qn(o->buffer[i]);
uint64_t t = Get_tn(o->buffer[i]);
int q_hit = 0, t_hit = 0;
k = h64_get(h, q);
q_hit = (k < kh_end(h) && kh_val(h, k) < m);
k = h64_get(h, t);
t_hit = (k < kh_end(h) && kh_val(h, k) < m);
if ((!q_hit && !t_hit) || (q_hit && t_hit)) continue;
if (!q_hit) {
k = h64_put(h, q, &absent);
if (absent) kh_val(h, k) = m;
}
if (!t_hit) {
k = h64_put(h, t, &absent);
if (absent) kh_val(h, k) = m;
}
}
}
}
}
for (j = 0; j < rs->total_reads; ++j) {
for (l = 0; l < 2; ++l) {
const ma_hit_t_alloc *o = &ov[l][j];
for (i = 0; i < (int)o->length; ++i) {
uint64_t q = Get_qn(o->buffer[i]);
uint64_t t = Get_tn(o->buffer[i]);
int q_hit = 0, t_hit = 0;
q_hit = (h64_get(h, q) < kh_end(h));
t_hit = (h64_get(h, t) < kh_end(h));
if (!q_hit && !t_hit) continue;
fwrite(Get_NAME(*rs, q), 1, Get_NAME_LENGTH(*rs, q), fp);
fwrite("\t", 1, 1, fp);
fprintf(fp, "%lu\t", (unsigned long)Get_READ_LENGTH(*rs, q));
fprintf(fp, "%d\t", Get_qs(o->buffer[i]));
fprintf(fp, "%d\t", Get_qe(o->buffer[i]));
fputs(o->buffer[i].rev? "-\t" : "+\t", fp);
fwrite(Get_NAME(*rs, t), 1, Get_NAME_LENGTH(*rs, t), fp);
fwrite("\t", 1, 1, fp);
fprintf(fp, "%lu\t", (unsigned long)Get_READ_LENGTH(*rs, t));
fprintf(fp, "%d\t", Get_ts(o->buffer[i]));
fprintf(fp, "%d\t%d\t%d\t%d\n", Get_te(o->buffer[i]), o->buffer[i].ml, o->buffer[i].bl, !l);
}
}
}
h64_destroy(h);
}
void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o)
{
int i, n;
char **list;
list = gv_read_list(o, &n);
ha_extract_print(rs, n_rounds, n, list);
for (i = 0; i < n; ++i) free(list[i]);
free(list);
}
+11858
View File
File diff suppressed because it is too large Load Diff
+72
View File
@@ -0,0 +1,72 @@
#ifndef __HIC__
#define __HIC__
#include <stdint.h>
#include "Overlaps.h"
#define kdq_clear(q) ((q)->count = (q)->front = 0)
#define kv_malloc(v, s) ((v).n = 0, (v).m = (s), MALLOC((v).a, (s)))
#define RC_0 0
#define RC_1 1
#define RC_2 2
hc_edge* get_hc_edge(hc_links* link, uint64_t src, uint64_t dest, uint64_t dir);
void push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t* d);
void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link);
void hic_benchmark(ma_ug_t *ug, asg_t* read_g);
typedef struct {
double w;
uint32_t id, occ;
///uint32_t *bid, bid_n;
ma_utg_t *u;
uint64_t l_d, r_d;
}chain_hic_w_type;
typedef struct {
size_t n, m;
chain_hic_w_type* a;
uint32_t max_bub_id;
uint32_t *chain_idx, u_n;
}chain_hic_warp;
typedef struct {
long long g_occ, b_occ;
uint64_t id;
uint8_t del;
}chain_w_type;
typedef struct {
uint32_t* index, round_id, n_round;
ma_ug_t* ug;
kvec_t(uint32_t) list;
kvec_t(uint32_t) num;
kvec_t(uint64_t) pathLen;
kvec_t(uint64_t) b_s_idx;
uint64_t s_bub, f_bub, b_bub, b_end_bub, tangle_bub, cross_bub, mess_bub;
uint32_t check_het;
asg_t *b_g;
ma_ug_t* b_ug;
kvec_t(chain_w_type) chain_weight;
chain_hic_warp c_w;
} bubble_type;
#define P_het(B) ((B).num.n)
#define M_het(B) ((B).num.n + 1)
// #define IF_BUB(ID, B) ((B).index[(ID)] < (B).num.n)
// #define IF_HET(ID, B) ((B).index[(ID)] == (B).num.n)
// #define IF_HOM(ID, B) ((B).index[(ID)] > (B).num.n)
#define IF_BUB(ID, B) ((B).index[(ID)] < (B).f_bub+1)
#define IF_HET(ID, B) ((B).index[(ID)] == (B).f_bub+1)
#define IF_HOM(ID, B) ((B).index[(ID)] > (B).f_bub+1)
#define Get_bub_num(RECORD) ((RECORD).num.n-1)
void get_bubbles(bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t* sink, uint32_t** a, uint32_t* n, uint64_t* pathBase);
int load_hc_links(hc_links* link, const char *fn);
void write_hc_links(hc_links* link, const char *fn);
void destory_bubbles(bubble_type* bub);
void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link);
void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub);
uint32_t connect_bub_occ(bubble_type* bub, uint32_t root_id, uint32_t check_het);
void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, uint32_t check_het);
void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint32_t is_end);
void set_b_utg_weight_flag(bubble_type* bub, buf_t* b, uint32_t v, uint8_t* vis_flag, uint32_t flag, uint32_t* occ);
#endif
+31 -3
View File
@@ -1,4 +1,4 @@
.TH hifiasm 1 "12 Apr 2020" "hifiasm-0.5.0" "Bioinformatics tools"
.TH hifiasm 1 "19 July 2020" "hifiasm-0.9 (r289)" "Bioinformatics tools"
.SH NAME
.PP
@@ -111,6 +111,11 @@ assembly.
.BI -r \ INT
Rounds of haplotype-aware error corrections [3]. This option affects all outputs of hifiasm.
.TP
.BI --min-hist-cnt \ INT
When analyzing the k-mer spectrum, ignore counts below
.IR INT .
.SS Assembly options
.TP
@@ -192,6 +197,20 @@ and do the assembly directly and quickly.
This might be helpful when users want to get an optimized assembly by multiple rounds of experiments
with different parameters.
.TP
.BI --pri-range \ INT1[,INT2]
Min and max coverage cutoff of primary contigs.
Keep contigs with coverage in this range at p_ctg.gfa.
Inferred automatically in default.
If INT2 is not specified, it is set to infinity.
Set -1 to disable.
.TP
.BI --lowQ \ INT
Output contig regions with >=INT% inconsistency to the bed file
with suffix
.B lowQ.bed
[70]. Set 0 to disable.
.SS Trio-partition options
@@ -252,6 +271,17 @@ Similarity threshold for duplicate haplotigs that should be purged [0.75].
.BI -O \ FLOAT
Min number of overlapped reads for duplicate haplotigs that should be purged [1].
.TP
.BI --purge-cov \ INT
Coverage upper bound of Purge-dups, which is inferred automatically in default.
If the coverage of a contig is higher than this bound, don't apply Purge-dups.
.TP
.BI --high-het \ INT
Enable this mode for high heterozygosity sample, which will increase running time.
For ordinary samples, no need to enable this mode [experimental, not stable].
.SS Debugging options
.TP 10
@@ -322,5 +352,3 @@ For each graph, hifiasm also outputs a simplified version without sequences for
the ease of visualization. Hifiasm keeps corrected reads and overlaps in three
binary files such as it can regenerate assembly graphs from the binary files
without redoing error correction.
+7 -3
View File
@@ -1,4 +1,5 @@
#include <stdio.h>
#include <assert.h>
#include "htab.h"
static void ha_hist_line(int c, int x, int exceed, int64_t cnt)
@@ -11,16 +12,19 @@ static void ha_hist_line(int c, int x, int exceed, int64_t cnt)
fprintf(stderr, " %lld\n", (long long)cnt);
}
int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het)
int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het)
{
const int hist_max = 100;
int i, start, low_i, max_i, max2_i, max3_i;
int64_t max, max2, max3, min;
// find the low point from the left
// determine the start point
assert(n_cnt > start_cnt);
*peak_het = -1;
start = cnt[1] > 0? 1 : 2;
low_i = start;
// find the low point from the left
low_i = start > start_cnt? start : start_cnt;
for (i = low_i + 1; i < n_cnt; ++i)
if (cnt[i] > cnt[i-1]) break;
low_i = i - 1;
+425 -29
View File
@@ -35,6 +35,9 @@ const unsigned char seq_nt4_table[256] = { // translate ACGT to 0123
void *ha_flt_tab;
ha_pt_t *ha_idx;
void *ha_flt_tab_hp;
ha_pt_t *ha_idx_hp;
void *ha_ct_table;
/***************************
* Yak specific parameters *
@@ -46,6 +49,7 @@ typedef struct {
int32_t pre;
int32_t n_thread;
int64_t chunk_size;
int adaLen;
} yak_copt_t;
void yak_copt_init(yak_copt_t *o)
@@ -71,7 +75,7 @@ typedef struct {
int n_shift, n_hashes;
uint8_t *b;
} yak_bf_t;
///in most cases, n_shift = 25, n_hashes = 4
yak_bf_t *yak_bf_init(int n_shift, int n_hashes)
{
yak_bf_t *b;
@@ -126,10 +130,12 @@ typedef struct {
typedef struct {
int k, pre, n_hash, n_shift;
uint64_t tot;
uint64_t tot; ///number of distinct k-mers
ha_ct1_t *h;
} ha_ct_t;
///for 0-th counting, k = 51, pre = 12, n_hash = 4, n_shift = 37
///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0
static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
{
ha_ct_t *h;
@@ -138,12 +144,15 @@ static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
CALLOC(h, 1);
h->k = k, h->pre = pre;
CALLOC(h->h, 1<<h->pre);
///i<h->pre = 4096
///it seems there is a large hash table h, consisting 4096 small hash tables
for (i = 0; i < 1<<h->pre; ++i)
h->h[i].h = yak_ct_init();
///for 0-th counting, enter here; used for bloom filter
if (n_hash > 0 && n_shift > h->pre) {
h->n_hash = n_hash, h->n_shift = n_shift;
for (i = 0; i < 1<<h->pre; ++i)
h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash);
h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash); ///h->n_shift = 37, h->pre = 12, h->n_hash = 4
}
return h;
}
@@ -173,15 +182,23 @@ static int ha_ct_insert_list(ha_ct_t *h, int create_new, int n, const uint64_t *
int j, mask = (1<<h->pre) - 1, n_ins = 0;
ha_ct1_t *g;
if (n == 0) return 0;
///corresponding small hash index
g = &h->h[a[0]&mask];
for (j = 0; j < n; ++j) {
int ins = 1, absent;
///x is a 64-bit word, h->pre=12
///all elements at a have the same low 12 bits
///so low 12 bits are not useful
uint64_t x = a[j] >> h->pre;
khint_t k;
if ((a[j]&mask) != (a[0]&mask)) continue;
if (create_new) {
///for 0-th counting, g->b = NULL
if (g->b)
ins = (yak_bf_insert(g->b, x) == h->n_hash);
///for 0-th counting, g->b = NULL
///x = the high 52 bits of a[j] + low 12 bits 0
///the low 12 bits are used for counting
if (ins) {
k = yak_ct_put(g->h, x << YAK_COUNTER_BITS | (g->b? 1 : 0), &absent);
if (absent) ++n_ins;
@@ -219,6 +236,8 @@ static void worker_ct_hist(void *data, long i, int tid) // callback for kt_for()
++cnt[kh_key(g, k)&YAK_MAX_COUNT];
}
///YAK_N_COUNTS is also 4096
///used for calculating k-mer histogram
static void ha_ct_hist(const ha_ct_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread)
{
hist_aux_t a;
@@ -226,6 +245,7 @@ static void ha_ct_hist(const ha_ct_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread
a.h = h;
memset(cnt, 0, YAK_N_COUNTS * sizeof(uint64_t));
CALLOC(a.cnt, n_thread);
///start 4096 threads
kt_for(n_thread, worker_ct_hist, &a, 1<<h->pre);
for (i = 0; i < YAK_N_COUNTS; ++i) cnt[i] = 0;
for (j = 0; j < n_thread; ++j)
@@ -265,6 +285,7 @@ static void ha_ct_shrink(ha_ct_t *h, int min, int max, int n_thread)
int i;
shrink_aux_t a;
a.h = h, a.min = min, a.max = max;
///still start 4096 threads
kt_for(n_thread, worker_ct_shrink, &a, 1<<h->pre);
for (i = 0, h->tot = 0; i < 1<<h->pre; ++i)
h->tot += kh_size(h->h[i].h);
@@ -306,6 +327,7 @@ static void worker_pt_gen(void *data, long i, int tid) // callback for kt_for()
int absent;
khint_t l;
l = yak_pt_put(b->h, kh_key(g, k) >> a->ct->pre << YAK_COUNTER_BITS, &absent);
///this should be the start index of kh_key's corresponding pos at ha_idxpos_t* a
kh_val(b->h, l) = b->n;
b->n += kh_key(g, k) & YAK_MAX_COUNT;
}
@@ -412,8 +434,12 @@ typedef struct {
ha_mz1_t *b;
} ch_buf_t;
///p = 12
static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer
{
///assign k-mer to one of the 4096 bins
///using low 12 bits for assigning
///so all elements at b have the same low 12 bits
int pre = y & ((1<<p) - 1);
ch_buf_t *b = &buf[pre];
if (b->n == b->m) {
@@ -425,6 +451,7 @@ static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k
static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y)
{
///assign minimizer to one of 4096 bins by low 12 bits
int pre = y->x & ((1<<p) - 1);
ch_buf_t *b = &buf[pre];
if (b->n == b->m) {
@@ -434,13 +461,17 @@ static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y)
b->b[b->n++] = *y;
}
///buf is the read block, k is the k-mer length, p = 12, len is the read length, seq is the read
static void count_seq_buf(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf
{
int i, l;
uint64_t x[4], mask = (1ULL<<k) - 1, shift = k - 1;
for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) {
int c = seq_nt4_table[(uint8_t)seq[i]];
///c = 00, 01, 10, 11
if (c < 4) { // not an "N" base
///x[0] & x[1] are the forward k-mer
///x[2] & x[3] are the reverse complementary k-mer
x[0] = (x[0] << 1 | (c&1)) & mask;
x[1] = (x[1] << 1 | (c>>1)) & mask;
x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift;
@@ -483,12 +514,13 @@ KSEQ_INIT(gzFile, gzread)
#define HAF_RS_WRITE_SEQ 0x8
#define HAF_RS_READ 0x10
#define HAF_CREATE_NEW 0x20
#define HAF_SKIP_READ 0x40
typedef struct { // global data structure for kt_pipeline()
const yak_copt_t *opt;
const void *flt_tab;
int flag, create_new, is_store;
uint64_t n_seq;
uint64_t n_seq; ///number of total reads
kseq_t *ks;
UC_Read ucr;
ha_ct_t *ct;
@@ -499,7 +531,8 @@ typedef struct { // global data structure for kt_pipeline()
typedef struct { // data structure for each step in kt_pipeline()
pl_data_t *p;
uint64_t n_seq0;
uint64_t n_seq0; ///the start index of current buffer block at R_INF
///sum_len = total bases, nk = number of k-mers
int n_seq, m_seq, sum_len, nk;
int *len;
char **seq;
@@ -514,15 +547,17 @@ static void worker_for_insert(void *data, long i, int tid) // callback for kt_fo
ch_buf_t *b = &s->buf[i];
if (s->p->pt)
b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b);
else
else///for 0-th count, go into here
b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a);
}
static void worker_for_mz(void *data, long i, int tid)
{
st_data_t *s = (st_data_t*)data;
///get the corresponding minimzer vector of this read
ha_mz1_v *b = &s->mz_buf[tid];
s->mz_buf[tid].n = 0;
///s->p->opt->w = 51, s->p->opt->k
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);
s->mz[i].n = s->mz[i].m = b->n;
MALLOC(s->mz[i].a, b->n);
@@ -540,6 +575,11 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
s->n_seq0 = p->n_seq;
if (p->rs_in && (p->flag & HAF_RS_READ)) {
while (p->n_seq < p->rs_in->total_reads) {
if((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU)
{
++p->n_seq;
continue;
}
int l;
recover_UC_Read(&p->ucr, p->rs_in, p->n_seq);
l = p->ucr.length;
@@ -559,12 +599,15 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
}
} else {
while ((ret = kseq_read(p->ks)) >= 0) {
int l = p->ks->seq.l;
int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen);
if(l <= 0) continue;
if (p->n_seq >= 1<<28) {
fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);
exit(1);
}
if (p->rs_out) {
///for 0-th count, just insert read length to R_INF, instead of read
if (p->flag & HAF_RS_WRITE_LEN) {
assert(p->n_seq == p->rs_out->total_reads);
ha_insert_read_len(p->rs_out, l, p->ks->name.l);
@@ -572,23 +615,25 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
int i, n_N;
assert(l == (int)p->rs_out->read_length[p->n_seq]);
for (i = n_N = 0; i < l; ++i) // count number of ambiguous bases
if (seq_nt4_table[(uint8_t)p->ks->seq.s[i]] >= 4)
if (seq_nt4_table[(uint8_t)p->ks->seq.s[i+p->opt->adaLen]] >= 4)
++n_N;
ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s, l, &p->rs_out->N_site[p->n_seq], n_N);
ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N);
memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);
}
}
///for 0-th count, insert both seq and length to local block
if (s->n_seq == s->m_seq) {
s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
REALLOC(s->len, s->m_seq);
REALLOC(s->seq, s->m_seq);
}
MALLOC(s->seq[s->n_seq], l);
memcpy(s->seq[s->n_seq], p->ks->seq.s, l);
memcpy(s->seq[s->n_seq], p->ks->seq.s+p->opt->adaLen, l);
s->len[s->n_seq++] = l;
++p->n_seq;
s->sum_len += l;
s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
///p->opt->chunk_size is the block max size
if (s->sum_len >= p->opt->chunk_size)
break;
}
@@ -596,18 +641,24 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
if (s->sum_len == 0) free(s);
else return s;
} else if (step == 1) { // step 2: extract k-mers
///s is the block of reads
st_data_t *s = (st_data_t*)in;
///for 0-th counting, n_pre = 4096
int i, n_pre = 1<<p->opt->pre, m;
// allocate the k-mer buffer
CALLOC(s->buf, n_pre);
m = (int)(s->nk * 1.2 / n_pre) + 1;
//pre-allocate memory for each of 4096 buffer
for (i = 0; i < n_pre; ++i) {
s->buf[i].m = m;
///for 0-th counting, p->pt = NULL
if (p->pt) MALLOC(s->buf[i].b, m);
else MALLOC(s->buf[i].a, m);
}
// fill the buffer
///for 0-th counting, p->opt->w == 1
if (p->opt->w == 1) { // enumerate all k-mers
///scan all reads
for (i = 0; i < s->n_seq; ++i) {
if (p->opt->is_HPC)
count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
@@ -618,18 +669,21 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
} else { // minimizers only
uint32_t j;
// compute minimizers
// s->n_seq is how many reads at this buffer
// s->mz && s->mz_buf are lists of minimzer vectors
CALLOC(s->mz, s->n_seq);
CALLOC(s->mz_buf, 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)
free(s->mz_buf[i].a);
free(s->mz_buf);
// insert minimizers
if (p->pt) {
if (p->pt) {///insert whole minimizer
for (i = 0; i < s->n_seq; ++i)
for (j = 0; j < s->mz[i].n; ++j)
pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]);
} else {
} else {///just insert the hash key of minimizer
for (i = 0; i < s->n_seq; ++i)
for (j = 0; j < s->mz[i].n; ++j)
ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
@@ -640,6 +694,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
}
free(s->mz);
}
///just clean seq
free(s->seq); free(s->len);
s->seq = 0, s->len = 0;
return s;
@@ -647,7 +702,9 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
st_data_t *s = (st_data_t*)in;
int i, n = 1<<p->opt->pre;
uint64_t n_ins = 0;
///for 0-th counting, p->pt = NULL
kt_for(p->opt->n_thread, worker_for_insert, s, n);
///n_ins is number of distinct k-mers
for (i = 0; i < n; ++i) {
n_ins += s->buf[i].n_ins;
if (p->pt) free(s->buf[i].b);
@@ -666,8 +723,51 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
return 0;
}
void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs)
{
int ret;
uint32_t i, m, pass, unpass;
gzFile fp = 0;
kseq_t *ks = NULL;
UC_Read ucr;
init_UC_Read(&ucr);
for (i = m = pass = unpass = 0; i < (uint32_t)asm_opt->num_reads; ++i)
{
if ((fp = gzopen(asm_opt->read_file_names[i], "r")) == 0) continue;
ks = kseq_init(fp);
while ((ret = kseq_read(ks)) >= 0)
{
int l = ks->seq.l;
if((l - asm_opt->adapterLen*2) <= 0) continue;
recover_UC_Read(&ucr, rs, m);
fprintf(stderr, "l: %d, ucr.length: %lld, asm_opt->adapterLen: %d\n",
l, ucr.length, asm_opt->adapterLen);
if(memcmp(ucr.seq, ks->seq.s+asm_opt->adapterLen, ucr.length) == 0)
{
pass++;
}
else
{
unpass++;
}
m++;
}
kseq_destroy(ks);
gzclose(fp);
ks = NULL;
fp = 0;
}
destory_UC_Read(&ucr);
fprintf(stderr, "[M::%s::# reads: %u, # pass: %u, # unpass: %u\n]", __func__, m, pass, unpass);
exit(1);
}
static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt_t *p0, ha_ct_t *c0, const void *flt_tab, All_reads *rs, int64_t *n_seq)
{
///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN|HAF_CREATE_NEW
int read_rs = (rs && (flag & HAF_RS_READ));
pl_data_t pl;
gzFile fp = 0;
@@ -676,23 +776,29 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
if (read_rs) {
pl.rs_in = rs;
init_UC_Read(&pl.ucr);
} else {
} else {///for 0-th counting, go into here
if ((fp = gzopen(fn, "r")) == 0) return 0;
pl.ks = kseq_init(fp);
}
///for 0-th counting, read all reads into pl.rs_out
if (rs && (flag & (HAF_RS_WRITE_LEN|HAF_RS_WRITE_SEQ)))
pl.rs_out = rs;
///for 0-th counting, flt_tab = NULL
///for 1-th counting, flt_tab = NULL
pl.flt_tab = flt_tab;
pl.opt = opt;
pl.flag = flag;
if (p0) {
if (p0) {///for 1-th counting, p0 = NULL
pl.pt = p0, pl.create_new = 0; // never create new elements in a position table
assert(p0->k == opt->k && p0->pre == opt->pre);
} else if (c0) {
pl.ct = c0, pl.create_new = !!(flag&HAF_CREATE_NEW);
assert(c0->k == opt->k && c0->pre == opt->pre);
} else {
} else {///for ft-th counting and 1-th counting, go into here
pl.create_new = 1; // alware create new elements if the count table is empty
///for 0-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 37
///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0
///building a large hash table consisting of 4096 small hash tables
pl.ct = ha_ct_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift);
}
kt_pipeline(3, worker_count, &pl, 3);
@@ -713,6 +819,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
yak_copt_t opt;
ha_ct_t *h = 0;
assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both
///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN
if (rs) {
if (flag & HAF_RS_WRITE_LEN)
init_All_reads(rs);
@@ -721,10 +828,16 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
}
yak_copt_init(&opt);
opt.k = asm_opt->k_mer_length;
///always 0
opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC);
///for ft-counting, shoud be 1
opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win;
///for ft-counting, shoud be 37
///for ha_pt_gen, shoud be 0
opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift;
opt.n_thread = asm_opt->thread_num;
opt.adaLen = asm_opt->adapterLen;
///asm_opt->num_reads is the number of fastq files
for (i = 0; i < asm_opt->num_reads; ++i)
h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, &n_seq);
if (h && opt.bf_shift > 0)
@@ -771,23 +884,63 @@ void ha_ft_destroy(void *h)
if (h) yak_ft_destroy((yak_ft_t*)h);
}
void debug_ct_index(void* q_ct_idx, void* r_ct_idx)
{
ha_ct_t* ct_idx = (ha_ct_t*)q_ct_idx;
yak_ct_t *g = NULL;
uint64_t i;
khint_t k;
for (i = 0; (int)i < 1<<ct_idx->pre; i++)
{
g = ct_idx->h[i].h;
for (k = 0; k < kh_end(g); ++k)
{
if (kh_exist(g, k))
{
int c = kh_key(g, k) & YAK_MAX_COUNT;
uint64_t hash = ((kh_key(g, k) >> ct_idx->pre)<<ct_idx->pre) | i;
int q = query_ct_index(r_ct_idx, hash);
if(q!=c)
{
fprintf(stderr, "ERROR:c: %d, q: %d\n", c, q);
}
}
}
}
}
/*************************
* High-level interfaces *
*************************/
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode)
{
yak_ft_t *flt_tab;
int64_t cnt[YAK_N_COUNTS];
int peak_hom, peak_het, cutoff;
int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0;
if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ;
ha_ct_t *h;
h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN, NULL, NULL, rs);
ha_ct_hist(h, cnt, asm_opt->thread_num);
peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
if (hom_cov) *hom_cov = peak_hom;
if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
cutoff = (int)(peak_hom * asm_opt->high_factor);
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs);
if((asm_opt->flag & HA_F_VERBOSE_GFA))
{
write_ct_index((void*)h, asm_opt->output_file_name);
// load_ct_index(&ha_ct_table, asm_opt->output_file_name);
// debug_ct_index((void*)h, ha_ct_table);
// debug_ct_index(ha_ct_table, (void*)h);
// ha_ct_destroy((ha_ct_t *)ha_ct_table);
}
if(!(ex_flag & HAF_SKIP_READ))
{
ha_ct_hist(h, cnt, asm_opt->thread_num);
peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, cnt, &peak_het);
if (hom_cov) *hom_cov = peak_hom;
if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
///in default, asm_opt->high_factor = 5.0
cutoff = (int)(peak_hom * asm_opt->high_factor);
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
}
ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num);
flt_tab = gen_hh(h);
ha_ct_destroy(h);
@@ -796,35 +949,42 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
return (void*)flt_tab;
}
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov)
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov)
{
int64_t cnt[YAK_N_COUNTS], tot_cnt;
int peak_hom, peak_het, i, extra_flag1, extra_flag2;
ha_ct_t *ct;
ha_pt_t *pt;
if (read_from_store) {
if (read_from_store) {///if reads have already been read
extra_flag1 = extra_flag2 = HAF_RS_READ;
} else if (rs->total_reads == 0) {
} else if (rs->total_reads == 0) {///if reads & length have not been scanned
extra_flag1 = HAF_RS_WRITE_LEN;
extra_flag2 = HAF_RS_WRITE_SEQ;
} else {
} else {///if length has been loaded but reads have not
extra_flag1 = HAF_RS_WRITE_SEQ;
extra_flag2 = HAF_RS_READ;
}
if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ;
ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
yak_realtime(), yak_cpu_usage(), (long)ct->tot);
ha_ct_hist(ct, cnt, asm_opt->thread_num);
fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]);
peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, cnt, &peak_het);
if (hom_cov) *hom_cov = peak_hom;
if (het_cov) *het_cov = peak_het;
if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
///here ha_ct_shrink is mostly used to remove k-mer appearing only 1 time
if (flt_tab == 0) {
int cutoff = (int)(peak_hom * asm_opt->high_factor);
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
if((extra_flag1 & HAF_SKIP_READ) && (extra_flag2 & HAF_SKIP_READ)) cutoff = YAK_MAX_COUNT - 1;
ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num);
for (i = 2, tot_cnt = 0; i <= cutoff; ++i) tot_cnt += cnt[i] * i;
} else {
///Note: here is just to remove minimizer appearing YAK_MAX_COUNT times
///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen
ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num);
for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i;
}
@@ -836,3 +996,239 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos);
return pt;
}
int query_ct_index(void* ct_idx, uint64_t hash)
{
ha_ct1_t *g = &(((ha_ct_t*)ct_idx)->h[hash & ((1ULL<<((ha_ct_t*)ct_idx)->pre) - 1)]);
khint_t k;
k = yak_ct_get(g->h, hash);
if (k == kh_end(g->h)) return 0;
return kh_key(g->h, k)&YAK_MAX_COUNT;
}
int write_ct_index(void *i_ct_idx, char* file_name)
{
char* gfa_name = (char*)malloc(strlen(file_name)+25);
sprintf(gfa_name, "%s.ct_flt", file_name);
FILE* fp = fopen(gfa_name, "w");
if (!fp) {
free(gfa_name);
return 0;
}
ha_ct_t* ct_idx = (ha_ct_t*)i_ct_idx;
int i;
ha_ct1_t *g;
fwrite(&ct_idx->k, sizeof(ct_idx->k), 1, fp);
fwrite(&ct_idx->pre, sizeof(ct_idx->pre), 1, fp);
fwrite(&ct_idx->n_hash, sizeof(ct_idx->n_hash), 1, fp);
fwrite(&ct_idx->n_shift, sizeof(ct_idx->n_shift), 1, fp);
fwrite(&ct_idx->tot, sizeof(ct_idx->tot), 1, fp);
for (i = 0; i < 1<<ct_idx->pre; i++)
{
g = &(ct_idx->h[i]);
yak_ct_save(g->h, fp);
}
fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
free(gfa_name);
fclose(fp);
return 1;
}
int load_ct_index(void **i_ct_idx, char* file_name)
{
char* gfa_name = (char*)malloc(strlen(file_name)+25);
sprintf(gfa_name, "%s.ct_flt", file_name);
FILE* fp = fopen(gfa_name, "r");
if (!fp) {
free(gfa_name);
return 0;
}
ha_ct_t** ct_idx = (ha_ct_t**)i_ct_idx;
double index_time = 0;
uint64_t flag = 0;
int i;
ha_ct_t *h = 0;
ha_ct1_t *g;
CALLOC(h, 1);
flag += fread(&h->k, sizeof(h->k), 1, fp);
flag += fread(&h->pre, sizeof(h->pre), 1, fp);
flag += fread(&h->n_hash, sizeof(h->n_hash), 1, fp);
flag += fread(&h->n_shift, sizeof(h->n_shift), 1, fp);
flag += fread(&h->tot, sizeof(h->tot), 1, fp);
CALLOC(h->h, 1<<h->pre);
index_time = yak_realtime();
for (i = 0; i < 1<<h->pre; ++i)
{
g = &(h->h[i]);
yak_ct_load(&(g->h), fp);
}
(*ct_idx) = h;
fprintf(stderr, "[M::%s::%.3f] ==> Loaded count table\n", __func__, yak_realtime() - index_time);
fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__);
free(gfa_name);
return 1;
}
int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
{
char* gfa_name = (char*)malloc(strlen(file_name)+25);
sprintf(gfa_name, "%s.pt_flt", file_name);
FILE* fp = fopen(gfa_name, "w");
if (!fp) {
free(gfa_name);
return 0;
}
yak_ft_t *ha_flt_tab = (yak_ft_t*)flt_tab;
if(ha_flt_tab)
{
fwrite("f", 1, 1, fp);
yak_ft_save(ha_flt_tab, fp);
}
if(ha_idx)
{
int i;
ha_pt1_t *g;
fwrite("h", 1, 1, fp);
fwrite(&ha_idx->k, sizeof(ha_idx->k), 1, fp);
fwrite(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp);
fwrite(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp);
fwrite(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp);
for (i = 0; i < 1<<ha_idx->pre; ++i)
{
g = &(ha_idx->h[i]);
yak_pt_save(g->h, fp);
fwrite(&g->n, sizeof(g->n), 1, fp);
fwrite(g->a, sizeof(ha_idxpos_t), g->n, fp);
}
}
fwrite(&opt->number_of_round, sizeof(opt->number_of_round), 1, fp);
fwrite(&opt->hom_cov, sizeof(opt->hom_cov), 1, fp);
fwrite(&opt->het_cov, sizeof(opt->het_cov), 1, fp);
fwrite(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp);
write_All_reads(r, gfa_name);
fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
free(gfa_name);
fclose(fp);
return 1;
}
int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
{
char* gfa_name = (char*)malloc(strlen(file_name)+25);
sprintf(gfa_name, "%s.pt_flt", file_name);
FILE* fp = fopen(gfa_name, "r");
if (!fp) {
free(gfa_name);
return 0;
}
ha_pt_t *ha_idx = NULL;
char mode = 0;
int f_flag, absent, i;
double index_time, index_s_time, pos_time, pos_s_time;
f_flag += fread(&mode, 1, 1, fp);
if(mode == 'f')
{
index_time = yak_realtime();
yak_ft_load((yak_ft_t **)r_flt_tab, fp);
f_flag += fread(&mode, 1, 1, fp);
fprintf(stderr, "[M::%s::%.3f] ==> Loaded flt table\n", __func__, yak_realtime()-index_time);
}
///insert using multiple threads???
if(mode == 'h')
{
pos_time = index_time = 0;
CALLOC(ha_idx, 1);
ha_pt1_t *g;
f_flag += fread(&ha_idx->k, sizeof(ha_idx->k), 1, fp);
f_flag += fread(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp);
f_flag += fread(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp);
f_flag += fread(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp);
CALLOC(ha_idx->h, 1<<ha_idx->pre);
for (i = 0; i < 1<<ha_idx->pre; ++i)
{
index_s_time = yak_realtime();
g = &(ha_idx->h[i]);
yak_pt_load(&(g->h), fp);
index_time += yak_realtime() - index_s_time;
pos_s_time = yak_realtime();
f_flag += fread(&g->n, sizeof(g->n), 1, fp);
MALLOC(g->a, g->n);
f_flag += fread(g->a, sizeof(ha_idxpos_t), g->n, fp);
pos_time += yak_realtime() - pos_s_time;
}
(*r_ha_idx) = ha_idx;
fprintf(stderr, "[M::%s::%.3f(index)/%.3f(pos)] ==> Loaded pos table\n", __func__, index_time, pos_time);
}
if(mode != 'h' && mode != 'f')
{
free(gfa_name);
fclose(fp);
return 0;
}
f_flag += fread(&absent, sizeof(absent), 1, fp);
if(absent != opt->number_of_round)
{
fprintf(stderr, "ERROR: different number of rounds!\n");
exit(1);
}
f_flag += fread(&opt->hom_cov, sizeof(opt->hom_cov), 1, fp);
f_flag += fread(&opt->het_cov, sizeof(opt->het_cov), 1, fp);
f_flag += fread(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp);
fclose(fp);
if(!load_All_reads(r, gfa_name))
{
free(gfa_name);
return 0;
}
memset(r->trio_flag, AMBIGU, r->total_reads*sizeof(uint8_t));
r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
for (i = 0; i < (long long)r->total_reads; i++)
{
init_ma_hit_t_alloc(&(r->paf[i]));
init_ma_hit_t_alloc(&(r->reverse_paf[i]));
}
fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__);
free(gfa_name);
return 1;
}
+19 -4
View File
@@ -6,7 +6,9 @@
#include "CommandLines.h"
typedef struct {
uint64_t x;
uint64_t x; ///x is the hash key
///rid is the read id, pos is the end pos of this minimizer, rev is the direction
///span is the length of this k-mer. For non-HPC k-mer, span may not be equal to k
uint64_t rid:28, pos:27, rev:1, span:8;
} ha_mz1_t;
@@ -25,15 +27,25 @@ typedef struct ha_abuf_s ha_abuf_t;
extern const unsigned char seq_nt4_table[256];
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
extern void *ha_flt_tab_hp;
extern ha_pt_t *ha_idx_hp;
extern void *ha_ct_table;
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov);
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode);
int ha_ft_isflt(const void *hh, uint64_t y);
void ha_ft_destroy(void *h);
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov);
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov);
void ha_pt_destroy(ha_pt_t *h);
const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);
int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
int write_ct_index(void *ct_idx, char* file_name);
int load_ct_index(void **ct_idx, char* file_name);
int query_ct_index(void* ct_idx, uint64_t hash);
ha_abuf_t *ha_abuf_init(void);
void ha_abuf_destroy(ha_abuf_t *ab);
uint64_t ha_abuf_mem(const ha_abuf_t *ab);
@@ -48,7 +60,9 @@ 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 ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het);
void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct);
int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het);
void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs);
static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function
{
@@ -76,6 +90,7 @@ static inline uint64_t yak_hash64_64(uint64_t key)
static inline uint64_t yak_hash_long(uint64_t x[4])
{
///compare forward k-mer and reverse complementary strand
int j = x[1] < x[3]? 0 : 1;
return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]);
}
+33
View File
@@ -133,6 +133,30 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
h->count = 0; \
} \
}
#define __KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \
SCOPE khint_t prefix##_save(HType *h, FILE* fp) { \
if (!h) return 0; \
khint_t n_buckets = (h->keys? 1U<<h->bits : 0U); \
fwrite(&n_buckets, sizeof(n_buckets), 1, fp); \
fwrite(&h->bits, sizeof(h->bits), 1, fp); \
fwrite(&h->count, sizeof(h->count), 1, fp); \
fwrite(h->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
fwrite(h->keys, sizeof(khkey_t), n_buckets, fp); \
return 1; \
} \
SCOPE khint_t prefix##_load(HType **h, FILE* fp) { \
(*h) = prefix##_init(); \
khint_t n_buckets; \
uint64_t flag = 0;\
flag += fread(&n_buckets, sizeof(n_buckets), 1, fp); \
flag += fread(&(*h)->bits, sizeof((*h)->bits), 1, fp); \
flag += fread(&(*h)->count, sizeof((*h)->count), 1, fp); \
(*h)->used = (khint32_t*)kmalloc(__kh_fsize(n_buckets) * sizeof(khint32_t)); \
(*h)->keys = (khkey_t*)kmalloc(n_buckets * sizeof(khkey_t)); \
flag += fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
flag += fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \
return 1; \
} \
#define __KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
SCOPE khint_t prefix##_getp(const HType *h, const khkey_t *key) { \
@@ -245,6 +269,7 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
#define KHASHL_INIT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
__KHASHL_TYPE(HType, khkey_t) \
__KHASHL_IMPL_BASIC(SCOPE, HType, prefix) \
__KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \
__KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
__KHASHL_IMPL_RESIZE(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
__KHASHL_IMPL_PUT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
@@ -264,6 +289,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
KHASHL_INIT(KH_LOCAL, HType, prefix##_s, HType##_s_bucket_t, prefix##_s_hash, prefix##_s_eq) \
SCOPE HType *prefix##_init(void) { return prefix##_s_init(); } \
SCOPE void prefix##_destroy(HType *h) { prefix##_s_destroy(h); } \
SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_s_save(h, fp); } \
SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_s_load(h, fp); } \
SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_s_resize(h, new_n_buckets); } \
SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_s_bucket_t t; t.key = key; return prefix##_s_getp(h, &t); } \
SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_s_del(h, k); } \
@@ -276,6 +303,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
KHASHL_INIT(KH_LOCAL, HType, prefix##_m, HType##_m_bucket_t, prefix##_m_hash, prefix##_m_eq) \
SCOPE HType *prefix##_init(void) { return prefix##_m_init(); } \
SCOPE void prefix##_destroy(HType *h) { prefix##_m_destroy(h); } \
SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_m_save(h, fp); } \
SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_m_load(h, fp); } \
SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_m_resize(h, new_n_buckets); } \
SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_m_bucket_t t; t.key = key; return prefix##_m_getp(h, &t); } \
SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_m_del(h, k); } \
@@ -287,6 +316,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
KHASHL_INIT(KH_LOCAL, HType, prefix##_cs, HType##_cs_bucket_t, __kh_cached_hash, prefix##_cs_eq) \
SCOPE HType *prefix##_init(void) { return prefix##_cs_init(); } \
SCOPE void prefix##_destroy(HType *h) { prefix##_cs_destroy(h); } \
SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_cs_save(h, fp); } \
SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_cs_load(h, fp); } \
SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cs_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cs_getp(h, &t); } \
SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cs_del(h, k); } \
SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cs_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cs_putp(h, &t, absent); }
@@ -297,6 +328,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
KHASHL_INIT(KH_LOCAL, HType, prefix##_cm, HType##_cm_bucket_t, __kh_cached_hash, prefix##_cm_eq) \
SCOPE HType *prefix##_init(void) { return prefix##_cm_init(); } \
SCOPE void prefix##_destroy(HType *h) { prefix##_cm_destroy(h); } \
SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_cm_save(h, fp); } \
SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_cm_load(h, fp); } \
SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cm_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cm_getp(h, &t); } \
SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cm_del(h, k); } \
SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cm_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cm_putp(h, &t, absent); }
+177
View File
@@ -0,0 +1,177 @@
#ifndef KSW2_H_
#define KSW2_H_
#include <stdint.h>
#define KSW_NEG_INF -0x40000000
#define KSW_EZ_SCORE_ONLY 0x01 // don't record alignment path/cigar
#define KSW_EZ_RIGHT 0x02 // right-align gaps
#define KSW_EZ_GENERIC_SC 0x04 // without this flag: match/mismatch only; last symbol is a wildcard
#define KSW_EZ_APPROX_MAX 0x08 // approximate max; this is faster with sse
#define KSW_EZ_APPROX_DROP 0x10 // approximate Z-drop; faster with sse
#define KSW_EZ_EXTZ_ONLY 0x40 // only perform extension
#define KSW_EZ_REV_CIGAR 0x80 // reverse CIGAR in the output
#define KSW_EZ_SPLICE_FOR 0x100
#define KSW_EZ_SPLICE_REV 0x200
#define KSW_EZ_SPLICE_FLANK 0x400
#ifdef __cplusplus
extern "C" {
#endif
typedef struct {
uint32_t max:31, zdropped:1;
int max_q, max_t; // max extension coordinate
int mqe, mqe_t; // max score when reaching the end of query
int mte, mte_q; // max score when reaching the end of target
int score; // max score reaching both ends; may be KSW_NEG_INF
int m_cigar, n_cigar;
int reach_end;
uint32_t *cigar;
} ksw_extz_t;
/**
* NW-like extension
*
* @param km memory pool, when used with kalloc
* @param qlen query length
* @param query query sequence with 0 <= query[i] < m
* @param tlen target length
* @param target target sequence with 0 <= target[i] < m
* @param m number of residue types
* @param mat m*m scoring mattrix in one-dimension array
* @param gapo gap open penalty; a gap of length l cost "-(gapo+l*gape)"
* @param gape gap extension penalty
* @param w band width (<0 to disable)
* @param zdrop off-diagonal drop-off to stop extension (positive; <0 to disable)
* @param flag flag (see KSW_EZ_* macros)
* @param ez (out) scores and cigar
*/
void ksw_extz(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez);
void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int flag, ksw_extz_t *ez);
void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
void ksw_extf2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t mch, int8_t mis, int8_t e, int w, int xdrop, ksw_extz_t *ez);
/**
* Global alignment
*
* (first 10 parameters identical to ksw_extz_sse())
* @param m_cigar (modified) max CIGAR length; feed 0 if cigar==0
* @param n_cigar (out) number of CIGAR elements
* @param cigar (out) BAM-encoded CIGAR; caller need to deallocate with kfree(km, )
*
* @return score of the alignment
*/
int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, const int8_t *mat);
int ksw_ll_i16(void *q, int tlen, const uint8_t *target, int gapo, int gape, int *qe, int *te);
#ifdef __cplusplus
}
#endif
/************************************
*** Private macros and functions ***
************************************/
#ifdef HAVE_KALLOC
#include "kalloc.h"
#else
#include <stdlib.h>
#define kmalloc(km, size) malloc((size))
#define kcalloc(km, count, size) calloc((count), (size))
#define krealloc(km, ptr, size) realloc((ptr), (size))
#define kfree(km, ptr) free((ptr))
#endif
static inline uint32_t *ksw_push_cigar(void *km, int *n_cigar, int *m_cigar, uint32_t *cigar, uint32_t op, int len)
{
if (*n_cigar == 0 || op != (cigar[(*n_cigar) - 1]&0xf)) {
if (*n_cigar == *m_cigar) {
*m_cigar = *m_cigar? (*m_cigar)<<1 : 4;
cigar = (uint32_t*)krealloc(km, cigar, (*m_cigar) << 2);
}
cigar[(*n_cigar)++] = len<<4 | op;
} else cigar[(*n_cigar)-1] += len<<4;
return cigar;
}
// In the backtrack matrix, value p[] has the following structure:
// bit 0-2: which type gets the max - 0 for H, 1 for E, 2 for F, 3 for \tilde{E} and 4 for \tilde{F}
// bit 3/0x08: 1 if a continuation on the E state (bit 5/0x20 for a continuation on \tilde{E})
// bit 4/0x10: 1 if a continuation on the F state (bit 6/0x40 for a continuation on \tilde{F})
static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int min_intron_len, const uint8_t *p, const int *off, const int *off_end, int n_col, int i0, int j0,
int *m_cigar_, int *n_cigar_, uint32_t **cigar_)
{ // p[] - lower 3 bits: which type gets the max; bit
int n_cigar = 0, m_cigar = *m_cigar_, i = i0, j = j0, r, state = 0;
uint32_t *cigar = *cigar_, tmp;
while (i >= 0 && j >= 0) { // at the beginning of the loop, _state_ tells us which state to check
int force_state = -1;
if (is_rot) {
r = i + j;
if (i < off[r]) force_state = 2;
if (off_end && i > off_end[r]) force_state = 1;
tmp = force_state < 0? p[(size_t)r * n_col + i - off[r]] : 0;
} else {
if (j < off[i]) force_state = 2;
if (off_end && j > off_end[i]) force_state = 1;
tmp = force_state < 0? p[(size_t)i * n_col + j - off[i]] : 0;
}
if (state == 0) state = tmp & 7; // if requesting the H state, find state one maximizes it.
else if (!(tmp >> (state + 2) & 1)) state = 0; // if requesting other states, _state_ stays the same if it is a continuation; otherwise, set to H
if (state == 0) state = tmp & 7; // TODO: probably this line can be merged into the "else if" line right above; not 100% sure
if (force_state >= 0) state = force_state;
if (state == 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 0, 1), --i, --j; // match
else if (state == 1 || (state == 3 && min_intron_len <= 0)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion
else if (state == 3 && min_intron_len > 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron
else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, 1), --j; // insertion
}
if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, min_intron_len > 0 && i >= min_intron_len? 3 : 2, i + 1); // first deletion
if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, j + 1); // first insertion
if (!is_rev)
for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR
tmp = cigar[i], cigar[i] = cigar[n_cigar-1-i], cigar[n_cigar-1-i] = tmp;
*m_cigar_ = m_cigar, *n_cigar_ = n_cigar, *cigar_ = cigar;
}
static inline void ksw_reset_extz(ksw_extz_t *ez)
{
ez->max_q = ez->max_t = ez->mqe_t = ez->mte_q = -1;
ez->max = 0, ez->score = ez->mqe = ez->mte = KSW_NEG_INF;
ez->n_cigar = 0, ez->zdropped = 0, ez->reach_end = 0;
}
static inline int ksw_apply_zdrop(ksw_extz_t *ez, int is_rot, int32_t H, int a, int b, int zdrop, int8_t e)
{
int r, t;
if (is_rot) r = a, t = b;
else r = a + b, t = a;
if (H > (int32_t)ez->max) {
ez->max = H, ez->max_t = t, ez->max_q = r - t;
} else if (t >= ez->max_t && r - t >= ez->max_q) {
int tl = t - ez->max_t, ql = (r - t) - ez->max_q, l;
l = tl > ql? tl - ql : ql - tl;
if (zdrop >= 0 && ez->max - H > zdrop + l * e) {
ez->zdropped = 1;
return 1;
}
}
return 0;
}
#endif
+305
View File
@@ -0,0 +1,305 @@
#include <string.h>
#include <assert.h>
#include "ksw2.h"
#ifdef __SSE2__
#include <emmintrin.h>
#ifdef KSW_SSE2_ONLY
#undef __SSE4_1__
#endif
#ifdef __SSE4_1__
#include <smmintrin.h>
#endif
#ifdef KSW_CPU_DISPATCH
#ifdef __SSE4_1__
void ksw_extz2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
#else
void ksw_extz2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
#endif
#else
void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
#endif // ~KSW_CPU_DISPATCH
{
#define __dp_code_block1 \
z = _mm_add_epi8(_mm_load_si128(&s[t]), qe2_); \
xt1 = _mm_load_si128(&x[t]); /* xt1 <- x[r-1][t..t+15] */ \
tmp = _mm_srli_si128(xt1, 15); /* tmp <- x[r-1][t+15] */ \
xt1 = _mm_or_si128(_mm_slli_si128(xt1, 1), x1_); /* xt1 <- x[r-1][t-1..t+14] */ \
x1_ = tmp; \
vt1 = _mm_load_si128(&v[t]); /* vt1 <- v[r-1][t..t+15] */ \
tmp = _mm_srli_si128(vt1, 15); /* tmp <- v[r-1][t+15] */ \
vt1 = _mm_or_si128(_mm_slli_si128(vt1, 1), v1_); /* vt1 <- v[r-1][t-1..t+14] */ \
v1_ = tmp; \
a = _mm_add_epi8(xt1, vt1); /* a <- x[r-1][t-1..t+14] + v[r-1][t-1..t+14] */ \
ut = _mm_load_si128(&u[t]); /* ut <- u[t..t+15] */ \
b = _mm_add_epi8(_mm_load_si128(&y[t]), ut); /* b <- y[r-1][t..t+15] + u[r-1][t..t+15] */
#define __dp_code_block2 \
z = _mm_max_epu8(z, b); /* z = max(z, b); this works because both are non-negative */ \
z = _mm_min_epu8(z, max_sc_); \
_mm_store_si128(&u[t], _mm_sub_epi8(z, vt1)); /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \
_mm_store_si128(&v[t], _mm_sub_epi8(z, ut)); /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \
z = _mm_sub_epi8(z, q_); \
a = _mm_sub_epi8(a, z); \
b = _mm_sub_epi8(b, z);
int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc;
int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX);
int32_t *H = 0, H0 = 0, last_H0_t = 0;
uint8_t *qr, *sf, *mem, *mem2 = 0;
__m128i q_, qe2_, zero_, flag1_, flag2_, flag8_, flag16_, sc_mch_, sc_mis_, sc_N_, m1_, max_sc_;
__m128i *u, *v, *x, *y, *s, *p = 0;
ksw_reset_extz(ez);
if (m <= 0 || qlen <= 0 || tlen <= 0) return;
zero_ = _mm_set1_epi8(0);
q_ = _mm_set1_epi8(q);
qe2_ = _mm_set1_epi8((q + e) * 2);
flag1_ = _mm_set1_epi8(1);
flag2_ = _mm_set1_epi8(2);
flag8_ = _mm_set1_epi8(0x08);
flag16_ = _mm_set1_epi8(0x10);
sc_mch_ = _mm_set1_epi8(mat[0]);
sc_mis_ = _mm_set1_epi8(mat[1]);
sc_N_ = mat[m*m-1] == 0? _mm_set1_epi8(-e) : _mm_set1_epi8(mat[m*m-1]);
m1_ = _mm_set1_epi8(m - 1); // wildcard
max_sc_ = _mm_set1_epi8(mat[0] + (q + e) * 2);
if (w < 0) w = tlen > qlen? tlen : qlen;
wl = wr = w;
tlen_ = (tlen + 15) / 16;
n_col_ = qlen < tlen? qlen : tlen;
n_col_ = ((n_col_ < w + 1? n_col_ : w + 1) + 15) / 16 + 1;
qlen_ = (qlen + 15) / 16;
for (t = 1, max_sc = mat[0], min_sc = mat[1]; t < m * m; ++t) {
max_sc = max_sc > mat[t]? max_sc : mat[t];
min_sc = min_sc < mat[t]? min_sc : mat[t];
}
if (-min_sc > 2 * (q + e)) return; // otherwise, we won't see any mismatches
mem = (uint8_t*)kcalloc(km, tlen_ * 6 + qlen_ + 1, 16);
u = (__m128i*)(((size_t)mem + 15) >> 4 << 4); // 16-byte aligned
v = u + tlen_, x = v + tlen_, y = x + tlen_, s = y + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16;
if (!approx_max) {
H = (int32_t*)kmalloc(km, tlen_ * 16 * 4);
for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF;
}
if (with_cigar) {
mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 16);
p = (__m128i*)(((size_t)mem2 + 15) >> 4 << 4);
off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2);
off_end = off + qlen + tlen - 1;
}
for (t = 0; t < qlen; ++t) qr[t] = query[qlen - 1 - t];
memcpy(sf, target, tlen);
for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) {
int st = 0, en = tlen - 1, st0, en0, st_, en_;
int8_t x1, v1;
uint8_t *qrr = qr + (qlen - 1 - r), *u8 = (uint8_t*)u, *v8 = (uint8_t*)v;
__m128i x1_, v1_;
// find the boundaries
if (st < r - qlen + 1) st = r - qlen + 1;
if (en > r) en = r;
if (st < (r-wr+1)>>1) st = (r-wr+1)>>1; // take the ceil
if (en > (r+wl)>>1) en = (r+wl)>>1; // take the floor
if (st > en) {
ez->zdropped = 1;
break;
}
st0 = st, en0 = en;
st = st / 16 * 16, en = (en + 16) / 16 * 16 - 1;
// set boundary conditions
if (st > 0) {
if (st - 1 >= last_st && st - 1 <= last_en)
x1 = ((uint8_t*)x)[st - 1], v1 = v8[st - 1]; // (r-1,s-1) calculated in the last round
else x1 = v1 = 0; // not calculated; set to zeros
} else x1 = 0, v1 = r? q : 0;
if (en >= r) ((uint8_t*)y)[r] = 0, u8[r] = r? q : 0;
// loop fission: set scores first
if (!(flag & KSW_EZ_GENERIC_SC)) {
for (t = st0; t <= en0; t += 16) {
__m128i sq, st, tmp, mask;
sq = _mm_loadu_si128((__m128i*)&sf[t]);
st = _mm_loadu_si128((__m128i*)&qrr[t]);
mask = _mm_or_si128(_mm_cmpeq_epi8(sq, m1_), _mm_cmpeq_epi8(st, m1_));
tmp = _mm_cmpeq_epi8(sq, st);
#ifdef __SSE4_1__
tmp = _mm_blendv_epi8(sc_mis_, sc_mch_, tmp);
tmp = _mm_blendv_epi8(tmp, sc_N_, mask);
#else
tmp = _mm_or_si128(_mm_andnot_si128(tmp, sc_mis_), _mm_and_si128(tmp, sc_mch_));
tmp = _mm_or_si128(_mm_andnot_si128(mask, tmp), _mm_and_si128(mask, sc_N_));
#endif
_mm_storeu_si128((__m128i*)((uint8_t*)s + t), tmp);
}
} else {
for (t = st0; t <= en0; ++t)
((uint8_t*)s)[t] = mat[sf[t] * m + qrr[t]];
}
// core loop
x1_ = _mm_cvtsi32_si128(x1);
v1_ = _mm_cvtsi32_si128(v1);
st_ = st / 16, en_ = en / 16;
assert(en_ - st_ + 1 <= n_col_);
if (!with_cigar) { // score only
for (t = st_; t <= en_; ++t) {
__m128i z, a, b, xt1, vt1, ut, tmp;
__dp_code_block1;
#ifdef __SSE4_1__
z = _mm_max_epi8(z, a); // z = z > a? z : a (signed)
#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8()
z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0;
z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative
#endif
__dp_code_block2;
#ifdef __SSE4_1__
_mm_store_si128(&x[t], _mm_max_epi8(a, zero_));
_mm_store_si128(&y[t], _mm_max_epi8(b, zero_));
#else
tmp = _mm_cmpgt_epi8(a, zero_);
_mm_store_si128(&x[t], _mm_and_si128(a, tmp));
tmp = _mm_cmpgt_epi8(b, zero_);
_mm_store_si128(&y[t], _mm_and_si128(b, tmp));
#endif
}
} else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment
__m128i *pr = p + (size_t)r * n_col_ - st_;
off[r] = st, off_end[r] = en;
for (t = st_; t <= en_; ++t) {
__m128i d, z, a, b, xt1, vt1, ut, tmp;
__dp_code_block1;
d = _mm_and_si128(_mm_cmpgt_epi8(a, z), flag1_); // d = a > z? 1 : 0
#ifdef __SSE4_1__
z = _mm_max_epi8(z, a); // z = z > a? z : a (signed)
tmp = _mm_cmpgt_epi8(b, z);
d = _mm_blendv_epi8(d, flag2_, tmp); // d = b > z? 2 : d
#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8()
z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0;
z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative
tmp = _mm_cmpgt_epi8(b, z);
d = _mm_or_si128(_mm_andnot_si128(tmp, d), _mm_and_si128(tmp, flag2_)); // d = b > z? 2 : d; emulating blendv
#endif
__dp_code_block2;
tmp = _mm_cmpgt_epi8(a, zero_);
_mm_store_si128(&x[t], _mm_and_si128(tmp, a));
d = _mm_or_si128(d, _mm_and_si128(tmp, flag8_)); // d = a > 0? 0x08 : 0
tmp = _mm_cmpgt_epi8(b, zero_);
_mm_store_si128(&y[t], _mm_and_si128(tmp, b));
d = _mm_or_si128(d, _mm_and_si128(tmp, flag16_)); // d = b > 0? 0x10 : 0
_mm_store_si128(&pr[t], d);
}
} else { // gap right-alignment
__m128i *pr = p + (size_t)r * n_col_ - st_;
off[r] = st, off_end[r] = en;
for (t = st_; t <= en_; ++t) {
__m128i d, z, a, b, xt1, vt1, ut, tmp;
__dp_code_block1;
d = _mm_andnot_si128(_mm_cmpgt_epi8(z, a), flag1_); // d = z > a? 0 : 1
#ifdef __SSE4_1__
z = _mm_max_epi8(z, a); // z = z > a? z : a (signed)
tmp = _mm_cmpgt_epi8(z, b);
d = _mm_blendv_epi8(flag2_, d, tmp); // d = z > b? d : 2
#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8()
z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0;
z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative
tmp = _mm_cmpgt_epi8(z, b);
d = _mm_or_si128(_mm_andnot_si128(tmp, flag2_), _mm_and_si128(tmp, d)); // d = z > b? d : 2; emulating blendv
#endif
__dp_code_block2;
tmp = _mm_cmpgt_epi8(zero_, a);
_mm_store_si128(&x[t], _mm_andnot_si128(tmp, a));
d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag8_)); // d = 0 > a? 0 : 0x08
tmp = _mm_cmpgt_epi8(zero_, b);
_mm_store_si128(&y[t], _mm_andnot_si128(tmp, b));
d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag16_)); // d = 0 > b? 0 : 0x10
_mm_store_si128(&pr[t], d);
}
}
if (!approx_max) { // find the exact max with a 32-bit score array
int32_t max_H, max_t;
// compute H[], max_H and max_t
if (r > 0) {
int32_t HH[4], tt[4], en1 = st0 + (en0 - st0) / 4 * 4, i;
__m128i max_H_, max_t_, qe_;
max_H = H[en0] = en0 > 0? H[en0-1] + u8[en0] - qe : H[en0] + v8[en0] - qe; // special casing the last element
max_t = en0;
max_H_ = _mm_set1_epi32(max_H);
max_t_ = _mm_set1_epi32(max_t);
qe_ = _mm_set1_epi32(q + e);
for (t = st0; t < en1; t += 4) { // this implements: H[t]+=v8[t]-qe; if(H[t]>max_H) max_H=H[t],max_t=t;
__m128i H1, tmp, t_;
H1 = _mm_loadu_si128((__m128i*)&H[t]);
t_ = _mm_setr_epi32(v8[t], v8[t+1], v8[t+2], v8[t+3]);
H1 = _mm_add_epi32(H1, t_);
H1 = _mm_sub_epi32(H1, qe_);
_mm_storeu_si128((__m128i*)&H[t], H1);
t_ = _mm_set1_epi32(t);
tmp = _mm_cmpgt_epi32(H1, max_H_);
#ifdef __SSE4_1__
max_H_ = _mm_blendv_epi8(max_H_, H1, tmp);
max_t_ = _mm_blendv_epi8(max_t_, t_, tmp);
#else
max_H_ = _mm_or_si128(_mm_and_si128(tmp, H1), _mm_andnot_si128(tmp, max_H_));
max_t_ = _mm_or_si128(_mm_and_si128(tmp, t_), _mm_andnot_si128(tmp, max_t_));
#endif
}
_mm_storeu_si128((__m128i*)HH, max_H_);
_mm_storeu_si128((__m128i*)tt, max_t_);
for (i = 0; i < 4; ++i)
if (max_H < HH[i]) max_H = HH[i], max_t = tt[i] + i;
for (; t < en0; ++t) { // for the rest of values that haven't been computed with SSE
H[t] += (int32_t)v8[t] - qe;
if (H[t] > max_H)
max_H = H[t], max_t = t;
}
} else H[0] = v8[0] - qe - qe, max_H = H[0], max_t = 0; // special casing r==0
// update ez
if (en0 == tlen - 1 && H[en0] > ez->mte)
ez->mte = H[en0], ez->mte_q = r - en;
if (r - st0 == qlen - 1 && H[st0] > ez->mqe)
ez->mqe = H[st0], ez->mqe_t = st0;
if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e)) break;
if (r == qlen + tlen - 2 && en0 == tlen - 1)
ez->score = H[tlen - 1];
} else { // find approximate max; Z-drop might be inaccurate, too.
if (r > 0) {
if (last_H0_t >= st0 && last_H0_t <= en0 && last_H0_t + 1 >= st0 && last_H0_t + 1 <= en0) {
int32_t d0 = v8[last_H0_t] - qe;
int32_t d1 = u8[last_H0_t + 1] - qe;
if (d0 > d1) H0 += d0;
else H0 += d1, ++last_H0_t;
} else if (last_H0_t >= st0 && last_H0_t <= en0) {
H0 += v8[last_H0_t] - qe;
} else {
++last_H0_t, H0 += u8[last_H0_t] - qe;
}
if ((flag & KSW_EZ_APPROX_DROP) && ksw_apply_zdrop(ez, 1, H0, r, last_H0_t, zdrop, e)) break;
} else H0 = v8[0] - qe - qe, last_H0_t = 0;
if (r == qlen + tlen - 2 && en0 == tlen - 1)
ez->score = H0;
}
last_st = st, last_en = en;
//for (t = st0; t <= en0; ++t) printf("(%d,%d)\t(%d,%d,%d,%d)\t%d\n", r, t, ((int8_t*)u)[t], ((int8_t*)v)[t], ((int8_t*)x)[t], ((int8_t*)y)[t], H[t]); // for debugging
}
kfree(km, mem);
if (!approx_max) kfree(km, H);
if (with_cigar) { // backtrack
int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR);
if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY)) {
ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
} else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int)ez->max) {
ez->reach_end = 1;
ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->mqe_t, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
} else if (ez->max_t >= 0 && ez->max_q >= 0) {
ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
}
kfree(km, mem2); kfree(km, off);
}
}
#endif // __SSE2__
+153 -2
View File
@@ -37,7 +37,11 @@ static inline int tq_shift(tiny_queue_t *q)
* @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)
{
{ ///in default, w = 51, k = 51, is_hpc = 1
/**
uint64_t x;
uint64_t rid:28, pos:27, rev:1, span:8;
**/
static const ha_mz1_t dummy = { UINT64_MAX, 0, 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;
@@ -45,8 +49,10 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
tiny_queue_t tq;
assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
///sizeof(ha_mz1_t) = 16
memset(buf, 0xff, w * 16);
memset(&tq, 0, sizeof(tiny_queue_t));
///len/w is the evaluated minimizer numbers
kv_resize(ha_mz1_t, *p, p->n + len/w);
for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
@@ -65,7 +71,11 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
tq_push(&tq, skip_len);
kmer_span += skip_len;
if (tq.count > k) kmer_span -= tq_shift(&tq);
} else kmer_span = l + 1 < k? l + 1 : k;
} else kmer_span = l + 1 < k? l + 1 : k;
///kmer_span should be used for HPC k-mer
///so for non-HPC k-mer, kmer_span should be k in any case?
///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
@@ -80,6 +90,12 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
}
} 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)
@@ -87,15 +103,26 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
for (j = 0; j < buf_pos; ++j)
if (min.x == buf[j].x && 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 (min.x >= buf[j].x) 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 (min.x >= buf[j].x) min = buf[j], min_pos = 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 (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
@@ -108,3 +135,127 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
if (min.x != UINT64_MAX)
kv_push(ha_mz1_t, *p, min);
}
void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf,
kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct)
{ ///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;
if(dbg_ct != NULL) dbg_ct->a.n = 0;
static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 };
uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0}, filtered;
int i, j, l, buf_pos, min_pos, kmer_span = 0;
ha_mz1_t buf[256], min = dummy;
tiny_queue_t tq;
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));
///sizeof(ha_mz1_t) = 16
memset(buf, 0xff, w * 16);
memset(&tq, 0, sizeof(tiny_queue_t));
///len/w is the evaluated minimizer numbers
kv_resize(ha_mz1_t, *p, p->n + len/w);
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
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;
y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
filtered = 0;
if(hf != 0) filtered = ha_ft_isflt(hf, y);
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 (hf == 0 || ha_ft_isflt(hf, y) == 0)
if(filtered == 0)
info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
if(k_flag != NULL) k_flag->a.a[i]++;
if(k_flag != NULL && filtered == 1) 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 (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
for (j = 0; j < buf_pos; ++j)
if (min.x == buf[j].x && 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 (min.x >= buf[j].x) 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 (min.x >= buf[j].x) min = buf[j], min_pos = 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 (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
for (j = 0; j <= buf_pos; ++j)
if (min.x == buf[j].x && 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);
}