X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=padding.c;h=3ee8087d4d054768f0f7d317e24acd208ac18683;hb=99b55f586822a33d012cd58a0e029d682ca4d70b;hp=b7a254ff69df23fbc99aa05c0c8ec162baa95251;hpb=9b352ed762c83f3df15cb385e0355930f6362074;p=samtools.git diff --git a/padding.c b/padding.c index b7a254f..3ee8087 100644 --- a/padding.c +++ b/padding.c @@ -166,11 +166,11 @@ int bam_pad2unpad(samfile_t *in, samfile_t *out, faidx_t *fai) uint32_t *cigar = bam1_cigar(b); n2 = 0; if (b->core.pos == 0 && b->core.tid >= 0 && strcmp(bam1_qname(b), h->target_name[b->core.tid]) == 0) { - // fprintf(stderr, "[depad] Found embedded reference %s\n", bam1_qname(b)); + // fprintf(stderr, "[depad] Found embedded reference '%s'\n", bam1_qname(b)); r_tid = b->core.tid; unpad_seq(b, &r); if (h->target_len[r_tid] != r.l) { - fprintf(stderr, "[depad] ERROR: (Padded) length of '%s' is %i in BAM header, but %id in embedded reference\n", bam1_qname(b), h->target_len[r_tid], r.l); + fprintf(stderr, "[depad] ERROR: (Padded) length of '%s' is %d in BAM header, but %ld in embedded reference\n", bam1_qname(b), h->target_len[r_tid], r.l); return -1; } if (fai) { @@ -186,8 +186,8 @@ int bam_pad2unpad(samfile_t *in, samfile_t *out, faidx_t *fai) // Show gaps as ASCII 45 fprintf(stderr, "[depad] ERROR: Embedded sequence and reference FASTA don't match for %s base %i, '%c' vs '%c'\n", h->target_name[b->core.tid], i+1, - r.s[i] ? bam_nt16_rev_table[r.s[i]] : 45, - q.s[i] ? bam_nt16_rev_table[q.s[i]] : 45); + r.s[i] ? bam_nt16_rev_table[(int)r.s[i]] : 45, + q.s[i] ? bam_nt16_rev_table[(int)q.s[i]] : 45); return -1; } } @@ -202,6 +202,7 @@ int bam_pad2unpad(samfile_t *in, samfile_t *out, faidx_t *fai) return -1; } else if (b->core.tid == r_tid) { ; // good case, reference available + //fprintf(stderr, "[depad] Have ref '%s' for read '%s'\n", h->target_name[b->core.tid], bam1_qname(b)); } else if (fai) { 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]); @@ -270,8 +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]; } + /* 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) { @@ -330,7 +366,7 @@ bam_header_t * fix_header(bam_header_t *old, faidx_t *fai) if (strlen(header->text) < header->l_text) { //fprintf(stderr, "[depad] Reallocating header buffer\n"); assert (newtext == header->text); - newtext = calloc(strlen(header->text) + 1, 1); + newtext = malloc(strlen(header->text) + 1); strcpy(newtext, header->text); free(header->text); header->text = newtext; @@ -414,11 +450,11 @@ int main_pad2unpad(int argc, char *argv[]) depad_end: // close files, free and return - free(fn_list); free(fn_out); - if (h != in->header) bam_header_destroy(h); + if (fai) fai_destroy(fai); + if (h != in->header) bam_header_destroy(h); samclose(in); samclose(out); - fai_destroy(fai); + free(fn_list); free(fn_out); return ret; } @@ -430,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"); @@ -438,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");