diff --git a/bwamem.c b/bwamem.c index 9950097..6c59f59 100644 --- a/bwamem.c +++ b/bwamem.c @@ -712,7 +712,7 @@ mem_alnreg_v mem_align1(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t * } // This routine is only used for the API purpose -mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, uint8_t *query, const mem_alnreg_t *ar) +mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, int l_query, uint8_t *query, const mem_alnreg_t *ar) { mem_aln_t a; int qb = ar->qb, qe = ar->qe, NM, score, is_rev; @@ -722,6 +722,18 @@ mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t * bwa_fix_xref(opt->mat, opt->q, opt->r, opt->w, bns, pac, (uint8_t*)query, &qb, &qe, &rb, &re); a.cigar = bwa_gen_cigar(opt->mat, opt->q, opt->r, opt->w, bns->l_pac, pac, qe - qb, (uint8_t*)&query[qb], rb, re, &score, &a.n_cigar, &NM); a.NM = NM; + if (qb != 0 || qe != l_query) { // add clipping to CIGAR + int clip5, clip3; + clip5 = is_rev? l_query - qe : qb; + clip3 = is_rev? qb : l_query - qe; + a.cigar = realloc(a.cigar, 4 * (a.n_cigar + 2)); + if (clip5) { + memmove(a.cigar+1, a.cigar, a.n_cigar * 4); + a.cigar[0] = clip5<<4|3; + ++a.n_cigar; + } + if (clip3) a.cigar[a.n_cigar++] = clip3<<4|3; + } pos = bns_depos(bns, rb < bns->l_pac? rb : re - 1, &is_rev); a.is_rev = is_rev; a.rid = bns_pos2rid(bns, pos); diff --git a/bwamem.h b/bwamem.h index 9996f6b..c2f124c 100644 --- a/bwamem.h +++ b/bwamem.h @@ -122,7 +122,7 @@ extern "C" { */ mem_alnreg_v mem_align1(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bns, const uint8_t *pac, int l_seq, char *seq); - mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, uint8_t *query, const mem_alnreg_t *ar); + mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, int l_query, uint8_t *query, const mem_alnreg_t *ar); /** * Infer the insert size distribution from interleaved alignment regions diff --git a/example.c b/example.c index 2fefde4..c0fede6 100644 --- a/example.c +++ b/example.c @@ -32,7 +32,7 @@ int main(int argc, char *argv[]) for (i = 0; i < ar.n; ++i) { // traverse each hit mem_aln_t a; if (ar.a[i].secondary >= 0) continue; // skip secondary alignments - a = mem_reg2aln(opt, idx->bns, idx->pac, (uint8_t*)ks->seq.s, &ar.a[i]); // get forward-strand position and CIGAR + a = mem_reg2aln(opt, idx->bns, idx->pac, ks->seq.l, (uint8_t*)ks->seq.s, &ar.a[i]); // get forward-strand position and CIGAR printf("%s\t%c\t%s\t%d\t%d\t", ks->name.s, "+-"[a.is_rev], idx->bns->anns[a.rid].name, a.pos, a.mapq); for (k = 0; k < a.n_cigar; ++k) printf("%d%c", a.cigar[k]>>4, "MIDSH"[a.cigar[k]&0xf]); diff --git a/main.c b/main.c index 0009fc6..d301bb4 100644 --- a/main.c +++ b/main.c @@ -4,7 +4,7 @@ #include "utils.h" #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.6.2-r308-beta" +#define PACKAGE_VERSION "0.6.2-r309-beta" #endif static int usage()