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

updated APIs abit

This commit is contained in:
Heng Li 2012-04-08 00:55:52 -04:00
parent 3b5a9e5595
commit 9c486fa41e
2 changed files with 42 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
bwa.c
View File

@ -24,6 +24,7 @@ struct bwa_idx_t {
struct bwa_buf_t { struct bwa_buf_t {
int max_buf; int max_buf;
bwa_pestat_t pes;
gap_stack_t *stack; gap_stack_t *stack;
gap_opt_t *opt; gap_opt_t *opt;
int *diff_tab; int *diff_tab;
@ -152,16 +153,15 @@ static void compute_NM(const uint8_t *pac, uint64_t l_pac, uint8_t *seq, int64_t
} }
} }
bwa_aln_t bwa_sa2aln(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, uint64_t sa, int n_gaps) void bwa_sa2aln(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, uint64_t sa, int n_gaps, bwa_aln_t *aln)
{ {
extern bwtint_t bwa_sa2pos(const bntseq_t *bns, const bwt_t *bwt, bwtint_t sapos, int len, int *strand); 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); 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, i, n_gap, n_mm; int strand, seq_len, i, n_gap, n_mm;
uint64_t pos3, pac_pos; uint64_t pos3, pac_pos;
uint8_t *s[2]; uint8_t *s[2];
bwa_aln_t aln;
memset(&aln, 0, sizeof(bwa_aln_t)); memset(aln, 0, sizeof(bwa_aln_t));
seq_len = strlen(seq); seq_len = strlen(seq);
if (seq_len<<1 > buf->max_buf) { if (seq_len<<1 > buf->max_buf) {
buf->max_buf = seq_len<<1; buf->max_buf = seq_len<<1;
@ -174,38 +174,41 @@ bwa_aln_t bwa_sa2aln(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, uint
s[0][i] = s[1][i] = nst_nt4_table[(int)seq[i]]; s[0][i] = s[1][i] = nst_nt4_table[(int)seq[i]];
seq_reverse(seq_len, s[1], 1); seq_reverse(seq_len, s[1], 1);
pac_pos = bwa_sa2pos(idx->bns, idx->bwt, sa, seq_len, &strand); pac_pos = bwa_sa2pos(idx->bns, idx->bwt, sa, seq_len, &strand);
if (strand) aln.flag |= 16; if (strand) aln->flag |= 16;
if (n_gaps) { // only for gapped alignment if (n_gaps) { // only for gapped alignment
int n_cigar; int n_cigar;
bwa_cigar_t *cigar16; bwa_cigar_t *cigar16;
cigar16 = bwa_refine_gapped_core(idx->bns->l_pac, idx->pac, seq_len, s[strand], &pac_pos, strand? n_gaps : -n_gaps, &n_cigar, 1); cigar16 = bwa_refine_gapped_core(idx->bns->l_pac, idx->pac, seq_len, s[strand], &pac_pos, strand? n_gaps : -n_gaps, &n_cigar, 1);
aln.n_cigar = n_cigar; aln->n_cigar = n_cigar;
aln.cigar = malloc(n_cigar * 4); aln->cigar = malloc(n_cigar * 4);
for (i = 0, pos3 = pac_pos; i < n_cigar; ++i) { for (i = 0, pos3 = pac_pos; i < n_cigar; ++i) {
int op = cigar16[i]>>14; int op = cigar16[i]>>14;
int len = cigar16[i]&0x3fff; int len = cigar16[i]&0x3fff;
if (op == 3) op = 4; // the 16-bit CIGAR is different from the 32-bit CIGAR if (op == 3) op = 4; // the 16-bit CIGAR is different from the 32-bit CIGAR
aln.cigar[i] = len<<4 | op; aln->cigar[i] = len<<4 | op;
if (op == 0 || op == 2) pos3 += len; if (op == 0 || op == 2) pos3 += len;
} }
free(cigar16); free(cigar16);
} else { // ungapped } else { // ungapped
aln.n_cigar = 1; aln->n_cigar = 1;
aln.cigar = malloc(4); aln->cigar = malloc(4);
aln.cigar[0] = seq_len<<4 | 0; aln->cigar[0] = seq_len<<4 | 0;
pos3 = pac_pos + seq_len; pos3 = pac_pos + seq_len;
} }
aln.n_n = bns_cnt_ambi(idx->bns, pac_pos, pos3 - pac_pos, &aln.ref_id); aln->n_n = bns_cnt_ambi(idx->bns, pac_pos, pos3 - pac_pos, &aln->ref_id);
aln.offset = pac_pos - idx->bns->anns[aln.ref_id].offset; aln->offset = pac_pos - idx->bns->anns[aln->ref_id].offset;
if (pos3 - idx->bns->anns[aln.ref_id].offset > idx->bns->anns[aln.ref_id].len) // read mapped beyond the end of a sequence if (pos3 - idx->bns->anns[aln->ref_id].offset > idx->bns->anns[aln->ref_id].len) // read mapped beyond the end of a sequence
aln.flag |= 4; // read unmapped aln->flag |= 4; // read unmapped
compute_NM(idx->pac, idx->bns->l_pac, s[strand], pac_pos, aln.n_cigar, aln.cigar, &n_mm, &n_gap); compute_NM(idx->pac, idx->bns->l_pac, s[strand], pac_pos, aln->n_cigar, aln->cigar, &n_mm, &n_gap);
aln.n_mm = n_mm; aln->n_mm = n_mm;
aln.n_gap = n_gap; aln->n_gap = n_gap;
return aln;
} }
bwa_one_t *bwa_se2(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, int gen_cigar) /************************
* Single-end alignment *
************************/
bwa_one_t *bwa_se(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, int gen_cigar)
{ {
bwa_one_t *one; bwa_one_t *one;
int best, cnt, i, seq_len; int best, cnt, i, seq_len;
@ -245,18 +248,25 @@ bwa_one_t *bwa_se2(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, int ge
one->mapQ = one->mapQs; one->mapQ = one->mapQs;
// compute CIGAR on request // compute CIGAR on request
one->one.ref_id = -1; one->one.ref_id = -1;
if (gen_cigar) one->one = bwa_sa2aln(idx, buf, seq, one->sa, one->which->n_gapo + one->which->n_gape); if (gen_cigar) bwa_sa2aln(idx, buf, seq, one->sa, one->which->n_gapo + one->which->n_gape, &one->one);
return one; return one;
} }
bwa_one_t *bwa_se(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq)
{
return bwa_se2(idx, buf, seq, 1);
}
void bwa_one_destroy(bwa_one_t *one) void bwa_one_destroy(bwa_one_t *one)
{ {
free(one->sai.sai); free(one->sai.sai);
free(one->one.cigar); free(one->one.cigar);
free(one); free(one);
} }
/************************
* Paired-end alignment *
************************/
void bwa_pestat(bwa_buf_t *buf, int n, bwa_one_t **o[2])
{
}
void bwa_pe(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq[2], bwa_one_t *o[2])
{
}

9
bwa.h
View File

@ -55,6 +55,11 @@ typedef struct {
bwa_aln_t one; bwa_aln_t one;
} bwa_one_t; } bwa_one_t;
typedef struct {
double avg, std, ap_prior;
uint64_t low, high, high_bayesian;
} bwa_pestat_t;
#ifdef __cplusplus #ifdef __cplusplus
extern "C" { extern "C" {
#endif #endif
@ -89,9 +94,9 @@ extern "C" {
* *
* @return An alignment * @return An alignment
*/ */
bwa_aln_t bwa_sa2aln(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, uint64_t sa, int n_gaps); void bwa_sa2aln(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, uint64_t sa, int n_gaps, bwa_aln_t *aln);
bwa_one_t *bwa_se(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq); bwa_one_t *bwa_se(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, int gen_cigar);
void bwa_one_destroy(bwa_one_t *one); void bwa_one_destroy(bwa_one_t *one);