From: pd3 Date: Tue, 23 Apr 2013 08:52:21 +0000 (-0700) Subject: Merge pull request #43 from peterjc/depad_bins X-Git-Url: https://git.donarmstrong.com/?a=commitdiff_plain;h=99b55f586822a33d012cd58a0e029d682ca4d70b;hp=4dde8f521add6fbc60641a280e38136c7b9e6ae5;p=samtools.git Merge pull request #43 from peterjc/depad_bins Update BIN values in 'samtools depad' & check in 'samtools index' --- diff --git a/bam_index.c b/bam_index.c index f916e04..2fac746 100644 --- a/bam_index.c +++ b/bam_index.c @@ -154,7 +154,7 @@ bam_index_t *bam_index_core(bamFile fp) bam_header_t *h; int i, ret; bam_index_t *idx; - uint32_t last_bin, save_bin; + uint32_t last_bin, save_bin, recalculated_bin; int32_t last_coor, last_tid, save_tid; bam1_core_t *c; uint64_t save_off, last_off, n_mapped, n_unmapped, off_beg, off_end, n_no_coor; @@ -193,6 +193,15 @@ bam_index_t *bam_index_core(bamFile fp) bam1_qname(b), last_coor, c->pos, c->tid+1); return NULL; } + if (c->tid >= 0) { + recalculated_bin = bam_reg2bin(c->pos, bam_calend(c, bam1_cigar(b))); + if (c->bin != recalculated_bin) { + fprintf(stderr, "[bam_index_core] read '%s' mapped at POS %d to %d has BIN %d but should be %d\n", + bam1_qname(b), c->pos + 1, bam_calend(c, bam1_cigar(b)), c->bin, recalculated_bin); + fprintf(stderr, "[bam_index_core] Fix it by using BAM->SAM->BAM to force a recalculation of the BIN field\n"); + return NULL; + } + } if (c->tid >= 0 && !(c->flag & BAM_FUNMAP)) insert_offset2(&idx->index2[b->core.tid], b, last_off); if (c->bin != last_bin) { // then possibly write the binning index if (save_bin != 0xffffffffu) // save_bin==0xffffffffu only happens to the first record 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");