padding.c

   1 #include <string.h>
   2 #include <assert.h>
   3 #include <unistd.h>
   4 #include "kstring.h"
   5 #include "sam_header.h"
   6 #include "sam.h"
   7 #include "bam.h"
   8 #include "faidx.h"
   9
  10 static void replace_cigar(bam1_t *b, int n, uint32_t *cigar)
  11 {
  12         if (n != b->core.n_cigar) {
  13                 int o = b->core.l_qname + b->core.n_cigar * 4;
  14                 if (b->data_len + (n - b->core.n_cigar) * 4 > b->m_data) {
  15                         b->m_data = b->data_len + (n - b->core.n_cigar) * 4;
  16                         kroundup32(b->m_data);
  17                         b->data = (uint8_t*)realloc(b->data, b->m_data);
  18                 }
  19                 memmove(b->data + b->core.l_qname + n * 4, b->data + o, b->data_len - o);
  20                 memcpy(b->data + b->core.l_qname, cigar, n * 4);
  21                 b->data_len += (n - b->core.n_cigar) * 4;
  22                 b->core.n_cigar = n;
  23         } else memcpy(b->data + b->core.l_qname, cigar, n * 4);
  24 }
  25
  26 #define write_cigar(_c, _n, _m, _v) do { \
  27                 if (_n == _m) { \
  28                         _m = _m? _m<<1 : 4; \
  29                         _c = (uint32_t*)realloc(_c, _m * 4); \
  30                 } \
  31                 _c[_n++] = (_v); \
  32         } while (0)
  33
  34 static void unpad_seq(bam1_t *b, kstring_t *s)
  35 {
  36         int k, j, i;
  37         int length;
  38         uint32_t *cigar = bam1_cigar(b);
  39         uint8_t *seq = bam1_seq(b);
  40         // b->core.l_qseq gives length of the SEQ entry (including soft clips, S)
  41         // We need the padded length after alignment from the CIGAR (excluding
  42         // soft clips S, but including pads from CIGAR D operations)
  43         length = 0;
  44         for (k = 0; k < b->core.n_cigar; ++k) {
  45                 int op, ol;
  46                 op= bam_cigar_op(cigar[k]);
  47                 ol = bam_cigar_oplen(cigar[k]);
  48                 if (op == BAM_CMATCH || op == BAM_CEQUAL || op == BAM_CDIFF || op == BAM_CDEL)
  49                         length += ol;
  50         }
  51         ks_resize(s, length);
  52         for (k = 0, s->l = 0, j = 0; k < b->core.n_cigar; ++k) {
  53                 int op, ol;
  54                 op = bam_cigar_op(cigar[k]);
  55                 ol = bam_cigar_oplen(cigar[k]);
  56                 if (op == BAM_CMATCH || op == BAM_CEQUAL || op == BAM_CDIFF) {
  57                         for (i = 0; i < ol; ++i, ++j) s->s[s->l++] = bam1_seqi(seq, j);
  58                 } else if (op == BAM_CSOFT_CLIP) {
  59                         j += ol;
  60                 } else if (op == BAM_CHARD_CLIP) {
  61                         /* do nothing */
  62                 } else if (op == BAM_CDEL) {
  63                         for (i = 0; i < ol; ++i) s->s[s->l++] = 0;
  64                 } else {
  65                         fprintf(stderr, "[depad] ERROR: Didn't expect CIGAR op %c in read %s\n", BAM_CIGAR_STR[op], bam1_qname(b));
  66                         assert(-1);
  67                 }
  68         }
  69         assert(length == s->l);
  70 }
  71
  72 int load_unpadded_ref(faidx_t *fai, char *ref_name, int ref_len, kstring_t *seq)
  73 {
  74         char base;
  75         char *fai_ref = 0;
  76         int fai_ref_len = 0, k;
  77
  78         fai_ref = fai_fetch(fai, ref_name, &fai_ref_len);
  79         if (fai_ref_len != ref_len) {
  80                 fprintf(stderr, "[depad] ERROR: FASTA sequence %s length %i, expected %i\n", ref_name, fai_ref_len, ref_len);
  81                 free(fai_ref);
  82                 return -1;
  83         }
  84         ks_resize(seq, ref_len);
  85         seq->l = 0;
  86         for (k = 0; k < ref_len; ++k) {
  87                 base = fai_ref[k];
  88                 if (base == '-' || base == '*') {
  89                         // Map gaps to null to match unpad_seq function
  90                         seq->s[seq->l++] = 0;
  91                 } else {
  92                         int i = bam_nt16_table[(int)base];
  93                         if (i == 0 || i==16) { // Equals maps to 0, anything unexpected to 16
  94                                 fprintf(stderr, "[depad] ERROR: Invalid character %c (ASCII %i) in FASTA sequence %s\n", base, (int)base, ref_name);
  95                                 free(fai_ref);
  96                                 return -1;
  97                         }
  98                         seq->s[seq->l++] = i;
  99                 }
 100         }
 101         assert(ref_len == seq->l);
 102         free(fai_ref);
 103         return 0;
 104 }
 105
 106 inline int * update_posmap(int *posmap, kstring_t ref)
 107 {
 108         int i, k;
 109         posmap = realloc(posmap, ref.m * sizeof(int));
 110         for (i = k = 0; i < ref.l; ++i) {
 111                 posmap[i] = k;
 112                 if (ref.s[i]) ++k;
 113         }
 114         return posmap;
 115 }
 116
 117 int bam_pad2unpad(samfile_t *in, samfile_t *out, faidx_t *fai)
 118 {
 119         bam_header_t *h;
 120         bam1_t *b;
 121         kstring_t r, q;
 122         int r_tid = -1;
 123         uint32_t *cigar2 = 0;
 124         int ret = 0, n2 = 0, m2 = 0, *posmap = 0;
 125
 126         b = bam_init1();
 127         r.l = r.m = q.l = q.m = 0; r.s = q.s = 0;
 128         int read_ret;
 129         h = in->header;
 130         while ((read_ret = samread(in, b)) >= 0) { // read one alignment from `in'
 131                 uint32_t *cigar = bam1_cigar(b);
 132                 n2 = 0;
 133                 if (b->core.pos == 0 && b->core.tid >= 0 && strcmp(bam1_qname(b), h->target_name[b->core.tid]) == 0) {
 134                         // fprintf(stderr, "[depad] Found embedded reference %s\n", bam1_qname(b));
 135                         r_tid = b->core.tid;
 136                         unpad_seq(b, &r);
 137                         if (h->target_len[r_tid] != r.l) {
 138                                 fprintf(stderr, "[depad] ERROR: (Padded) length of %s is %i in BAM header, but %ld in embedded reference\n", bam1_qname(b), h->target_len[r_tid], r.l);
 139                                 return -1;
 140                         }
 141                         if (fai) {
 142                                 // Check the embedded reference matches the FASTA file
 143                                 if (load_unpadded_ref(fai, h->target_name[b->core.tid], h->target_len[b->core.tid], &q)) return -1;
 144                                 assert(r.l == q.l);
 145                                 int i;
 146                                 for (i = 0; i < r.l; ++i) {
 147                                         if (r.s[i] != q.s[i]) {
 148                                                 // Show gaps as ASCII 45
 149                                                 fprintf(stderr, "[depad] ERROR: Embedded sequence and reference FASTA don't match for %s base %i, '%c' vs '%c'\n",
 150                                                         h->target_name[b->core.tid], i+1,
 151                                                         r.s[i] ? bam_nt16_rev_table[r.s[i]] : 45,
 152                                                         q.s[i] ? bam_nt16_rev_table[q.s[i]] : 45);
 153                                                 return -1;
 154                                         }
 155                                 }
 156                         }
 157                         write_cigar(cigar2, n2, m2, bam_cigar_gen(b->core.l_qseq, BAM_CMATCH));
 158                         replace_cigar(b, n2, cigar2);
 159                         posmap = update_posmap(posmap, r);
 160                 } else if (b->core.n_cigar > 0) {
 161                         int i, k, op;
 162                         if (b->core.tid < 0) {
 163                                 fprintf(stderr, "[depad] ERROR: Read %s has CIGAR but no RNAME\n", bam1_qname(b));
 164                                 return -1;
 165                         } else if (b->core.tid == r_tid) {
 166                                 ; // good case, reference available
 167                         } else if (fai) {
 168                                 if (load_unpadded_ref(fai, h->target_name[b->core.tid], h->target_len[b->core.tid], &r)) return -1;
 169                                 posmap = update_posmap(posmap, r);
 170                                 r_tid = b->core.tid;
 171                                 // fprintf(stderr, "[depad] Loaded %s from FASTA file\n", h->target_name[b->core.tid]);
 172                         } else {
 173                                 fprintf(stderr, "[depad] ERROR: Missing %s embedded reference sequence (and no FASTA file)\n", h->target_name[b->core.tid]);
 174                                 return -1;
 175                         }
 176                         unpad_seq(b, &q);
 177                         if (bam_cigar_op(cigar[0]) == BAM_CSOFT_CLIP) {
 178                                 write_cigar(cigar2, n2, m2, cigar[0]);
 179                         } else if (bam_cigar_op(cigar[0]) == BAM_CHARD_CLIP) {
 180                                 write_cigar(cigar2, n2, m2, cigar[0]);
 181                                 if (b->core.n_cigar > 2 && bam_cigar_op(cigar[1]) == BAM_CSOFT_CLIP) {
 182                                         write_cigar(cigar2, n2, m2, cigar[1]);
 183                                 }
 184                         }
 185                         /* Determine CIGAR operator for each base in the aligned read */
 186                         for (i = 0, k = b->core.pos; i < q.l; ++i, ++k)
 187                                 q.s[i] = q.s[i]? (r.s[k]? BAM_CMATCH : BAM_CINS) : (r.s[k]? BAM_CDEL : BAM_CPAD);
 188                         /* Include any pads if starts with an insert */
 189                         if (q.s[0] == BAM_CINS) {
 190                                 for (k = 0; k+1 < b->core.pos && !r.s[b->core.pos - k - 1]; ++k);
 191                                 if (k) write_cigar(cigar2, n2, m2, bam_cigar_gen(k, BAM_CPAD));
 192                         }
 193                         /* Count consecutive CIGAR operators to turn into a CIGAR string */
 194                         for (i = k = 1, op = q.s[0]; i < q.l; ++i) {
 195                                 if (op != q.s[i]) {
 196                                         write_cigar(cigar2, n2, m2, bam_cigar_gen(k, op));
 197                                         op = q.s[i]; k = 1;
 198                                 } else ++k;
 199                         }
 200                         write_cigar(cigar2, n2, m2, bam_cigar_gen(k, op));
 201                         if (bam_cigar_op(cigar[b->core.n_cigar-1]) == BAM_CSOFT_CLIP) {
 202                                 write_cigar(cigar2, n2, m2, cigar[b->core.n_cigar-1]);
 203                         } else if (bam_cigar_op(cigar[b->core.n_cigar-1]) == BAM_CHARD_CLIP) {
 204                                 if (b->core.n_cigar > 2 && bam_cigar_op(cigar[b->core.n_cigar-2]) == BAM_CSOFT_CLIP) {
 205                                         write_cigar(cigar2, n2, m2, cigar[b->core.n_cigar-2]);
 206                                 }
 207                                 write_cigar(cigar2, n2, m2, cigar[b->core.n_cigar-1]);
 208                         }
 209                         /* Remove redundant P operators between M/X/=/D operators, e.g. 5M2P10M -> 15M */
 210                         int pre_op, post_op;
 211                         for (i = 2; i < n2; ++i)
 212                                 if (bam_cigar_op(cigar2[i-1]) == BAM_CPAD) {
 213                                         pre_op = bam_cigar_op(cigar2[i-2]);
 214                                         post_op = bam_cigar_op(cigar2[i]);
 215                                         /* Note don't need to check for X/= as code above will use M only */
 216                                         if ((pre_op == BAM_CMATCH || pre_op == BAM_CDEL) && (post_op == BAM_CMATCH || post_op == BAM_CDEL)) {
 217                                                 /* This is a redundant P operator */
 218                                                 cigar2[i-1] = 0; // i.e. 0M
 219                                                 /* If had same operator either side, combine them in post_op */
 220                                                 if (pre_op == post_op) {
 221                                                         /* If CIGAR M, could treat as simple integers since BAM_CMATCH is zero*/
 222                                                         cigar2[i] = bam_cigar_gen(bam_cigar_oplen(cigar2[i-2]) + bam_cigar_oplen(cigar2[i]), post_op);
 223                                                         cigar2[i-2] = 0; // i.e. 0M
 224                                                 }
 225                                         }
 226                                 }
 227                         /* Remove the zero'd operators (0M) */
 228                         for (i = k = 0; i < n2; ++i)
 229                                 if (cigar2[i]) cigar2[k++] = cigar2[i];
 230                         n2 = k;
 231                         replace_cigar(b, n2, cigar2);
 232                         b->core.pos = posmap[b->core.pos];
 233                 }
 234                 samwrite(out, b);
 235         }
 236         if (read_ret < -1) {
 237                 fprintf(stderr, "[depad] truncated file.\n");
 238                 ret = 1;
 239         }
 240         free(r.s); free(q.s); free(posmap);
 241         bam_destroy1(b);
 242         return ret;
 243 }
 244
 245 static int usage(int is_long_help);
 246
 247 int main_pad2unpad(int argc, char *argv[])
 248 {
 249         samfile_t *in = 0, *out = 0;
 250         faidx_t *fai = 0;
 251         int c, is_bamin = 1, compress_level = -1, is_bamout = 1, is_long_help = 0;
 252         char in_mode[5], out_mode[5], *fn_out = 0, *fn_list = 0, *fn_ref = 0;
 253         int ret=0;
 254
 255         /* parse command-line options */
 256         strcpy(in_mode, "r"); strcpy(out_mode, "w");
 257         while ((c = getopt(argc, argv, "Sso:u1T:?")) >= 0) {
 258                 switch (c) {
 259                 case 'S': is_bamin = 0; break;
 260                 case 's': assert(compress_level == -1); is_bamout = 0; break;
 261                 case 'o': fn_out = strdup(optarg); break;
 262                 case 'u': assert(is_bamout == 1); compress_level = 0; break;
 263                 case '1': assert(is_bamout == 1); compress_level = 1; break;
 264                 case 'T': fn_ref = strdup(optarg); is_bamin = 0; break;
 265                 case '?': is_long_help = 1; break;
 266                 default: return usage(is_long_help);
 267                 }
 268         }
 269         if (argc == optind) return usage(is_long_help);
 270
 271         if (is_bamin) strcat(in_mode, "b");
 272         if (is_bamout) strcat(out_mode, "b");
 273         strcat(out_mode, "h");
 274         if (compress_level >= 0) {
 275                 char tmp[2];
 276                 tmp[0] = compress_level + '0'; tmp[1] = '\0';
 277                 strcat(out_mode, tmp);
 278         }
 279
 280         // Load FASTA reference (also needed for SAM -> BAM if missing header)
 281         if (fn_ref) {
 282                 fn_list = samfaipath(fn_ref);
 283                 fai = fai_load(fn_ref);
 284         }
 285         // open file handlers
 286         if ((in = samopen(argv[optind], in_mode, fn_list)) == 0) {
 287                 fprintf(stderr, "[depad] failed to open \"%s\" for reading.\n", argv[optind]);
 288                 ret = 1;
 289                 goto depad_end;
 290         }
 291         if (in->header == 0) {
 292                 fprintf(stderr, "[depad] failed to read the header from \"%s\".\n", argv[optind]);
 293                 ret = 1;
 294                 goto depad_end;
 295         }
 296         /* TODO - The reference sequence lengths in the BAM + SAM headers should be updated */
 297         if ((out = samopen(fn_out? fn_out : "-", out_mode, in->header)) == 0) {
 298                 fprintf(stderr, "[depad] failed to open \"%s\" for writing.\n", fn_out? fn_out : "standard output");
 299                 ret = 1;
 300                 goto depad_end;
 301         }
 302
 303         // Do the depad
 304         ret = bam_pad2unpad(in, out, fai);
 305
 306 depad_end:
 307         // close files, free and return
 308         free(fn_list); free(fn_out);
 309         samclose(in);
 310         samclose(out);
 311         fai_destroy(fai);
 312         return ret;
 313 }
 314
 315 static int usage(int is_long_help)
 316 {
 317         fprintf(stderr, "\n");
 318         fprintf(stderr, "Usage:   samtools depad <in.bam>\n\n");
 319         fprintf(stderr, "Options: -s       output is SAM (default is BAM)\n");
 320         fprintf(stderr, "         -S       input is SAM (default is BAM)\n");
 321         fprintf(stderr, "         -u       uncompressed BAM output (can't use with -s)\n");
 322         fprintf(stderr, "         -1       fast compression BAM output (can't use with -s)\n");
 323         fprintf(stderr, "         -T FILE  reference sequence file [null]\n");
 324         fprintf(stderr, "         -o FILE  output file name [stdout]\n");
 325         fprintf(stderr, "         -?       longer help\n");
 326         fprintf(stderr, "\n");
 327         if (is_long_help)
 328                 fprintf(stderr, "Notes:\n\
 329 \n\
 330   1. Requires embedded reference sequences (before the reads for that reference),\n\
 331      with the future aim to also support a FASTA padded reference sequence file.\n\
 332 \n\
 333   2. The input padded alignment read's CIGAR strings must not use P or I operators.\n\
 334 \n");
 335         return 1;
 336 }