X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=padding.c;h=a4066122833620e8a68abec0368458e096048d8b;hb=507c3e9db23c24b97d15bd1b34f321b6d21ecaee;hp=033a7916237f94314bcf9661880a3fcde4604156;hpb=0ab78ddb98818366103c11330526b8a2ea33a757;p=samtools.git diff --git a/padding.c b/padding.c index 033a791..a406612 100644 --- a/padding.c +++ b/padding.c @@ -30,9 +30,21 @@ static void replace_cigar(bam1_t *b, int n, uint32_t *cigar) static void unpad_seq(bam1_t *b, kstring_t *s) { int k, j, i; + int length; uint32_t *cigar = bam1_cigar(b); uint8_t *seq = bam1_seq(b); - ks_resize(s, b->core.l_qseq); + // b->core.l_qseq gives length of the SEQ entry (including soft clips, S) + // We need the padded length after alignment from the CIGAR (excluding + // soft clips S, but including pads from CIGAR D operations) + length = 0; + for (k = 0; k < b->core.n_cigar; ++k) { + int op, ol; + op= bam_cigar_op(cigar[k]); + ol = bam_cigar_oplen(cigar[k]); + if (op == BAM_CMATCH || op == BAM_CEQUAL || op == BAM_CDIFF || op == BAM_CDEL) + length += ol; + } + ks_resize(s, length); for (k = 0, s->l = 0, j = 0; k < b->core.n_cigar; ++k) { int op, ol; op = bam_cigar_op(cigar[k]); @@ -43,13 +55,14 @@ static void unpad_seq(bam1_t *b, kstring_t *s) j += ol; } else if (op == BAM_CHARD_CLIP) { /* do nothing */ - } else if (op == BAM_CDEL || op == BAM_CPAD) { + } else if (op == BAM_CDEL) { for (i = 0; i < ol; ++i) s->s[s->l++] = 0; } else { fprintf(stderr, "[depad] ERROR: Didn't expect CIGAR op %c in read %s\n", BAM_CIGAR_STR[op], bam1_qname(b)); assert(-1); } } + assert(length == s->l); } int bam_pad2unpad(bamFile in, bamFile out) @@ -97,19 +110,22 @@ int bam_pad2unpad(bamFile in, bamFile out) return -1; } unpad_seq(b, &q); - if (bam_cigar_op(cigar[0]) == BAM_CSOFT_CLIP) write_cigar(cigar2, n2, m2, cigar[0]); - if (bam_cigar_op(cigar[0]) == BAM_CHARD_CLIP) { + if (bam_cigar_op(cigar[0]) == BAM_CSOFT_CLIP) { + write_cigar(cigar2, n2, m2, cigar[0]); + } else if (bam_cigar_op(cigar[0]) == BAM_CHARD_CLIP) { write_cigar(cigar2, n2, m2, cigar[0]); if (b->core.n_cigar > 2 && bam_cigar_op(cigar[1]) == BAM_CSOFT_CLIP) { write_cigar(cigar2, n2, m2, cigar[1]); } } - /* Include any pads if starts with an insert */ - for (k = 0; k+1 < b->core.pos && !r.s[b->core.pos - k - 1]; ++k); - if (k) write_cigar(cigar2, n2, m2, bam_cigar_gen(k, BAM_CPAD)); /* Determine CIGAR operator for each base in the aligned read */ for (i = 0, k = b->core.pos; i < q.l; ++i, ++k) q.s[i] = q.s[i]? (r.s[k]? BAM_CMATCH : BAM_CINS) : (r.s[k]? BAM_CDEL : BAM_CPAD); + /* Include any pads if starts with an insert */ + if (q.s[0] == BAM_CINS) { + for (k = 0; k+1 < b->core.pos && !r.s[b->core.pos - k - 1]; ++k); + if (k) write_cigar(cigar2, n2, m2, bam_cigar_gen(k, BAM_CPAD)); + } /* Count consecutive CIGAR operators to turn into a CIGAR string */ for (i = k = 1, op = q.s[0]; i < q.l; ++i) { if (op != q.s[i]) { @@ -118,17 +134,33 @@ int bam_pad2unpad(bamFile in, bamFile out) } else ++k; } write_cigar(cigar2, n2, m2, bam_cigar_gen(k, op)); - if (bam_cigar_op(cigar[b->core.n_cigar-1]) == BAM_CSOFT_CLIP) write_cigar(cigar2, n2, m2, cigar[b->core.n_cigar-1]); - if (bam_cigar_op(cigar[b->core.n_cigar-1]) == BAM_CHARD_CLIP) { + if (bam_cigar_op(cigar[b->core.n_cigar-1]) == BAM_CSOFT_CLIP) { + write_cigar(cigar2, n2, m2, cigar[b->core.n_cigar-1]); + } else if (bam_cigar_op(cigar[b->core.n_cigar-1]) == BAM_CHARD_CLIP) { if (b->core.n_cigar > 2 && bam_cigar_op(cigar[b->core.n_cigar-2]) == BAM_CSOFT_CLIP) { write_cigar(cigar2, n2, m2, cigar[b->core.n_cigar-2]); } write_cigar(cigar2, n2, m2, cigar[b->core.n_cigar-1]); } - /* Remove redundant P operators between M operators, e.g. 5M2P10M -> 15M */ + /* Remove redundant P operators between M/X/=/D operators, e.g. 5M2P10M -> 15M */ + int pre_op, post_op; for (i = 2; i < n2; ++i) - if (bam_cigar_op(cigar2[i]) == BAM_CMATCH && bam_cigar_op(cigar2[i-1]) == BAM_CPAD && bam_cigar_op(cigar2[i-2]) == BAM_CMATCH) - cigar2[i] += cigar2[i-2], cigar2[i-2] = cigar2[i-1] = 0; + if (bam_cigar_op(cigar2[i-1]) == BAM_CPAD) { + pre_op = bam_cigar_op(cigar2[i-2]); + post_op = bam_cigar_op(cigar2[i]); + /* Note don't need to check for X/= as code above will use M only */ + if ((pre_op == BAM_CMATCH || pre_op == BAM_CDEL) && (post_op == BAM_CMATCH || post_op == BAM_CDEL)) { + /* This is a redundant P operator */ + cigar2[i-1] = 0; // i.e. 0M + /* If had same operator either side, combine them in post_op */ + if (pre_op == post_op) { + /* If CIGAR M, could treat as simple integers since BAM_CMATCH is zero*/ + cigar2[i] = bam_cigar_gen(bam_cigar_oplen(cigar2[i-2]) + bam_cigar_oplen(cigar2[i]), post_op); + cigar2[i-2] = 0; // i.e. 0M + } + } + } + /* Remove the zero'd operators (0M) */ for (i = k = 0; i < n2; ++i) if (cigar2[i]) cigar2[k++] = cigar2[i]; n2 = k; @@ -143,13 +175,15 @@ int bam_pad2unpad(bamFile in, bamFile out) return 0; } +static int usage(int is_long_help); + int main_pad2unpad(int argc, char *argv[]) { bamFile in, out; + int is_long_help = 1; int result=0; if (argc == 1) { - fprintf(stderr, "Usage: samtools depad \n"); - return 1; + return usage(is_long_help); } in = strcmp(argv[1], "-")? bam_open(argv[1], "r") : bam_dopen(fileno(stdin), "r"); out = bam_dopen(fileno(stdout), "w"); @@ -157,3 +191,29 @@ int main_pad2unpad(int argc, char *argv[]) bam_close(in); bam_close(out); return result; } + +static int usage(int is_long_help) +{ + fprintf(stderr, "\n"); + fprintf(stderr, "Usage: samtools depad \n\n"); + fprintf(stderr, "Currently there are no optional arguments.\n"); + //TODO - These are the arguments I think make sense to support: + //fprintf(stderr, "Usage: samtools depad [options] |\n\n"); + //fprintf(stderr, "Options: -b output BAM\n"); + //fprintf(stderr, " -S input is SAM\n"); + //fprintf(stderr, " -u uncompressed BAM output (force -b)\n"); + //fprintf(stderr, " -1 fast compression (force -b)\n"); + //fprintf(stderr, " -@ INT number of BAM compression threads [0]\n"); + //fprintf(stderr, " -T FILE reference sequence file (force -S) [null]\n"); + //fprintf(stderr, " -o FILE output file name [stdout]\n"); + //fprintf(stderr, " -? longer help\n"); + fprintf(stderr, "\n"); + if (is_long_help) + fprintf(stderr, "Notes:\n\ +\n\ + 1. Requires embedded reference sequences (before the reads for that reference).\n\ +\n\ + 2. The input padded alignment read's CIGAR strings must not use P or I operators.\n\ +\n"); + return 1; +}