diff --git a/bwa.c b/bwa.c index 4e2775f..8e99f18 100644 --- a/bwa.c +++ b/bwa.c @@ -24,6 +24,7 @@ struct bwa_idx_t { struct bwa_buf_t { int max_buf; + bwa_pestat_t pes; gap_stack_t *stack; gap_opt_t *opt; 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 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; uint64_t pos3, pac_pos; 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); if (seq_len<<1 > buf->max_buf) { 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]]; seq_reverse(seq_len, s[1], 1); 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 int n_cigar; 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); - aln.n_cigar = n_cigar; - aln.cigar = malloc(n_cigar * 4); + aln->n_cigar = n_cigar; + aln->cigar = malloc(n_cigar * 4); for (i = 0, pos3 = pac_pos; i < n_cigar; ++i) { int op = cigar16[i]>>14; int len = cigar16[i]&0x3fff; 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; } free(cigar16); } else { // ungapped - aln.n_cigar = 1; - aln.cigar = malloc(4); - aln.cigar[0] = seq_len<<4 | 0; + aln->n_cigar = 1; + aln->cigar = malloc(4); + aln->cigar[0] = seq_len<<4 | 0; pos3 = pac_pos + seq_len; } - 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; - 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 - 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_gap = n_gap; - return aln; + 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; + 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 + 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_gap = n_gap; } -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; 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; // compute CIGAR on request 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; } -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) { free(one->sai.sai); free(one->one.cigar); 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]) +{ +} diff --git a/bwa.h b/bwa.h index 5a75c13..e8172da 100644 --- a/bwa.h +++ b/bwa.h @@ -55,6 +55,11 @@ typedef struct { bwa_aln_t one; } bwa_one_t; +typedef struct { + double avg, std, ap_prior; + uint64_t low, high, high_bayesian; +} bwa_pestat_t; + #ifdef __cplusplus extern "C" { #endif @@ -89,9 +94,9 @@ extern "C" { * * @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);