From 66154ff5d2021bd7b1610bd26552399df9f212e5 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sat, 7 Apr 2012 01:25:39 -0400 Subject: [PATCH] towards refined gapped; unfinished --- bwa.c | 22 ++++++++++++++++++++++ bwa.h | 7 +++++++ bwase.c | 14 +++++++------- 3 files changed, 36 insertions(+), 7 deletions(-) diff --git a/bwa.c b/bwa.c index ab0afea..3c6fbfd 100644 --- a/bwa.c +++ b/bwa.c @@ -124,3 +124,25 @@ bwa_alnpre_t *bwa_aln_pre(const bwa_idx_t *idx, bwa_aux_t *aux, const char *seq, s[i] = s[i] > 3? 4 : 3 - s[i]; return (bwa_alnpre_t*)bwt_match_gap(idx->bwt, seq_len, s, w, seq_len <= aux->opt->seed_len? 0 : seed_w, &opt2, n_aln, aux->stack); } +/* +bwa_aln_t bwa_sa2aln(const bwa_idx_t *idx, bwa_aux_t *aux, const char *seq, uint64_t sa, int n_gaps) +{ + extern bwtint_t bwa_sa2pos(const bntseq_t *bns, const bwt_t *bwt, bwtint_t sapos, int len, int *strand); + extern bwa_cigar_t *bwa_refine_gapped_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const uint8_t *seq, bwtint_t *_pos, int ext, int *n_cigar, int is_end_correct); + int strand, seq_len, n_cigar; + uint64_t pos; + uint8_t *s[2]; + bwa_aln_t aln; + bwa_cigar_t *cigar16; + + seq_len = strlen(seq); + if (seq_len<<1 > aux->max_buf) { + aux->max_buf = seq_len<<1; + kroundup32(aux->max_buf); + aux->buf = realloc(aux->buf, aux->max_buf); + } + pos = bwa_sa2pos(idx->bns, idx->bwt, sa, seq_len, &strand); + cigar16 = bwa_refine_gapped_core(idx->bns->l_pac, idx->pac, len, s[strand], &pos, strand? n_gaps : -n_gaps, &n_cigar, 1); + return aln; +} +*/ diff --git a/bwa.h b/bwa.h index a69f9d2..01eddcf 100644 --- a/bwa.h +++ b/bwa.h @@ -25,6 +25,13 @@ typedef struct { uint64_t k, l; } bwa_alnpre_t; +typedef struct { + uint32_t n_cigar:15, gap:8, mm:8, strand:1; + uint32_t ref_id; + uint64_t offset; + uint32_t *cigar; +} bwa_aln_t; + extern bwa_opt_t bwa_def_opt; #ifdef __cplusplus diff --git a/bwase.c b/bwase.c index e2754cd..c88cae5 100644 --- a/bwase.c +++ b/bwase.c @@ -110,8 +110,8 @@ bwtint_t bwa_sa2pos(const bntseq_t *bns, const bwt_t *bwt, bwtint_t sapos, int l pos_f = bns_depos(bns, bwt_sa(bwt, sapos), &is_rev); // pos_f *strand = !is_rev; /* NB: For gapped alignment, pacpos may not be correct, which will be fixed - * in refine_gapped_core(). This line also determines the way "x" is - * calculated in refine_gapped_core() when (ext < 0 && is_end == 0). */ + * in bwa_refine_gapped_core(). This line also determines the way "x" is + * calculated in bwa_refine_gapped_core() when (ext < 0 && is_end == 0). */ if (is_rev) pos_f = pos_f + 1 < len? 0 : pos_f - len + 1; // mapped to the forward strand return pos_f; // FIXME: it is possible that pos_f < bns->anns[ref_id].offset } @@ -160,7 +160,7 @@ void bwa_cal_pac_pos(const bntseq_t *bns, const char *prefix, int n_seqs, bwa_se * forward strand. This happens when p->pos is calculated by * bwa_cal_pac_pos(). is_end_correct==0 if (*pos) gives the correct * coordinate. This happens only for color-converted alignment. */ -static bwa_cigar_t *refine_gapped_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const ubyte_t *seq, bwtint_t *_pos, +bwa_cigar_t *bwa_refine_gapped_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const ubyte_t *seq, bwtint_t *_pos, int ext, int *n_cigar, int is_end_correct) { bwa_cigar_t *cigar = 0; @@ -320,12 +320,12 @@ void bwa_refine_gapped(const bntseq_t *bns, int n_seqs, bwa_seq_t *seqs, ubyte_t bwt_multi1_t *q = s->multi + j; int n_cigar; if (q->gap == 0) continue; - q->cigar = refine_gapped_core(bns->l_pac, pacseq, s->len, q->strand? s->rseq : s->seq, &q->pos, + q->cigar = bwa_refine_gapped_core(bns->l_pac, pacseq, s->len, q->strand? s->rseq : s->seq, &q->pos, (q->strand? 1 : -1) * q->gap, &n_cigar, 1); q->n_cigar = n_cigar; } if (s->type == BWA_TYPE_NO_MATCH || s->type == BWA_TYPE_MATESW || s->n_gapo == 0) continue; - s->cigar = refine_gapped_core(bns->l_pac, pacseq, s->len, s->strand? s->rseq : s->seq, &s->pos, + s->cigar = bwa_refine_gapped_core(bns->l_pac, pacseq, s->len, s->strand? s->rseq : s->seq, &s->pos, (s->strand? 1 : -1) * (s->n_gapo + s->n_gape), &s->n_cigar, 1); } @@ -338,13 +338,13 @@ void bwa_refine_gapped(const bntseq_t *bns, int n_seqs, bwa_seq_t *seqs, ubyte_t int n_cigar; if (q->gap == 0) continue; free(q->cigar); - q->cigar = refine_gapped_core(bns->l_pac, ntpac, s->len, q->strand? s->rseq : s->seq, &q->pos, + q->cigar = bwa_refine_gapped_core(bns->l_pac, ntpac, s->len, q->strand? s->rseq : s->seq, &q->pos, (q->strand? 1 : -1) * q->gap, &n_cigar, 0); q->n_cigar = n_cigar; } if (s->type != BWA_TYPE_NO_MATCH && s->cigar) { // update cigar again free(s->cigar); - s->cigar = refine_gapped_core(bns->l_pac, ntpac, s->len, s->strand? s->rseq : s->seq, &s->pos, + s->cigar = bwa_refine_gapped_core(bns->l_pac, ntpac, s->len, s->strand? s->rseq : s->seq, &s->pos, (s->strand? 1 : -1) * (s->n_gapo + s->n_gape), &s->n_cigar, 0); } }