zzh
/
bwa
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

towards refined gapped; unfinished

This commit is contained in:
Heng Li 2012-04-07 01:25:39 -04:00
parent a1abfe9977
commit 66154ff5d2
3 changed files with 36 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
bwa.c
View File

@ -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;
}
*/

7
bwa.h
View File

@ -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

14
bwase.c
View File

@ -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);
}
}