X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=padding.c;h=3ee8087d4d054768f0f7d317e24acd208ac18683;hb=28391e5898804ce6b805016d8c676fdf61442eb3;hp=a8da562d0872e8b0930bb2bdd3dcc9abb58bab7b;hpb=f5c347706ae6ecce3e1794cb3141981a3c45397e;p=samtools.git diff --git a/padding.c b/padding.c index a8da562..3ee8087 100644 --- a/padding.c +++ b/padding.c @@ -271,40 +271,43 @@ int bam_pad2unpad(samfile_t *in, samfile_t *out, faidx_t *fai) if (cigar2[i]) cigar2[k++] = cigar2[i]; n2 = k; replace_cigar(b, n2, cigar2); - b->core.pos = posmap[b->core.pos]; - if (b->core.mtid < 0 || b->core.mpos < 0) { - /* Nice case, no mate to worry about*/ - // fprintf(stderr, "[depad] Read '%s' mate not mapped\n", bam1_qname(b)); - /* TODO - Warning if FLAG says mate should be mapped? */ - /* Clean up funny input where mate position is given but mate reference is missing: */ - b->core.mtid = -1; - b->core.mpos = -1; - } else if (b->core.mtid == b->core.tid) { - /* Nice case, same reference */ - // fprintf(stderr, "[depad] Read '%s' mate mapped to same ref\n", bam1_qname(b)); - b->core.mpos = posmap[b->core.mpos]; - } else { - /* Nasty case, Must load alternative posmap */ - // fprintf(stderr, "[depad] Loading reference '%s' temporarily\n", h->target_name[b->core.mtid]); - if (!fai) { - fprintf(stderr, "[depad] ERROR: Needed reference %s sequence for mate (and no FASTA file)\n", h->target_name[b->core.mtid]); - return -1; - } - /* Temporarily load the other reference sequence */ - if (load_unpadded_ref(fai, h->target_name[b->core.mtid], h->target_len[b->core.mtid], &r)) { - fprintf(stderr, "[depad] ERROR: Failed to load '%s' from reference FASTA\n", h->target_name[b->core.mtid]); - return -1; - } - posmap = update_posmap(posmap, r); - b->core.mpos = posmap[b->core.mpos]; - /* Restore the reference and posmap*/ - if (load_unpadded_ref(fai, h->target_name[b->core.tid], h->target_len[b->core.tid], &r)) { - fprintf(stderr, "[depad] ERROR: Failed to load '%s' from reference FASTA\n", h->target_name[b->core.tid]); - return -1; - } - posmap = update_posmap(posmap, r); + } + /* Even unmapped reads can have a POS value, e.g. if their mate was mapped */ + if (b->core.pos != -1) b->core.pos = posmap[b->core.pos]; + if (b->core.mtid < 0 || b->core.mpos < 0) { + /* Nice case, no mate to worry about*/ + // fprintf(stderr, "[depad] Read '%s' mate not mapped\n", bam1_qname(b)); + /* TODO - Warning if FLAG says mate should be mapped? */ + /* Clean up funny input where mate position is given but mate reference is missing: */ + b->core.mtid = -1; + b->core.mpos = -1; + } else if (b->core.mtid == b->core.tid) { + /* Nice case, same reference */ + // fprintf(stderr, "[depad] Read '%s' mate mapped to same ref\n", bam1_qname(b)); + b->core.mpos = posmap[b->core.mpos]; + } else { + /* Nasty case, Must load alternative posmap */ + // fprintf(stderr, "[depad] Loading reference '%s' temporarily\n", h->target_name[b->core.mtid]); + if (!fai) { + fprintf(stderr, "[depad] ERROR: Needed reference %s sequence for mate (and no FASTA file)\n", h->target_name[b->core.mtid]); + return -1; + } + /* Temporarily load the other reference sequence */ + if (load_unpadded_ref(fai, h->target_name[b->core.mtid], h->target_len[b->core.mtid], &r)) { + fprintf(stderr, "[depad] ERROR: Failed to load '%s' from reference FASTA\n", h->target_name[b->core.mtid]); + return -1; } + posmap = update_posmap(posmap, r); + b->core.mpos = posmap[b->core.mpos]; + /* Restore the reference and posmap*/ + if (load_unpadded_ref(fai, h->target_name[b->core.tid], h->target_len[b->core.tid], &r)) { + fprintf(stderr, "[depad] ERROR: Failed to load '%s' from reference FASTA\n", h->target_name[b->core.tid]); + return -1; + } + posmap = update_posmap(posmap, r); } + /* Most reads will have been moved so safest to always recalculate the BIN value */ + b->core.bin = bam_reg2bin(b->core.pos, bam_calend(&b->core, bam1_cigar(b))); samwrite(out, b); } if (read_ret < -1) { @@ -463,7 +466,7 @@ static int usage(int is_long_help) fprintf(stderr, " -S input is SAM (default is BAM)\n"); fprintf(stderr, " -u uncompressed BAM output (can't use with -s)\n"); fprintf(stderr, " -1 fast compression BAM output (can't use with -s)\n"); - fprintf(stderr, " -T FILE reference sequence file [null]\n"); + fprintf(stderr, " -T FILE padded reference sequence file [null]\n"); fprintf(stderr, " -o FILE output file name [stdout]\n"); fprintf(stderr, " -? longer help\n"); fprintf(stderr, "\n"); @@ -471,7 +474,7 @@ static int usage(int is_long_help) fprintf(stderr, "Notes:\n\ \n\ 1. Requires embedded reference sequences (before the reads for that reference),\n\ - with the future aim to also support a FASTA padded reference sequence file.\n\ + or ideally a FASTA file of the padded reference sequences (via the -T argument).\n\ \n\ 2. The input padded alignment read's CIGAR strings must not use P or I operators.\n\ \n");