pysam/pysam_util.c

   1 #include <ctype.h>
   2 #include <assert.h>
   3 #include "bam.h"
   4 #include "khash.h"
   5 #include "ksort.h"
   6 #include "bam_endian.h"
   7 #include "knetfile.h"
   8 #include "pysam_util.h"
   9 #include "errmod.h" // for pysam_dump
  10
  11 #ifndef inline
  12 #define inline __inline
  13 #endif
  14
  15 // Definition of pysamerr
  16 #include "stdio.h"
  17 FILE * pysamerr = NULL;
  18
  19 FILE * pysam_set_stderr( FILE * f )
  20 {
  21   pysamerr = f;
  22   return f;
  23 }
  24
  25 // #######################################################
  26 // utility routines to avoid using callbacks in bam_fetch
  27 // taken from bam_index.c
  28 // The order of the following declarations is important.
  29 // #######################################################
  30
  31 typedef struct
  32 {
  33   uint64_t u, v;
  34 } pair64_t;
  35
  36 #define pair64_lt(a,b) ((a).u < (b).u)
  37
  38 typedef struct {
  39         uint32_t m, n;
  40         pair64_t *list;
  41 } bam_binlist_t;
  42
  43 typedef struct {
  44         int32_t n, m;
  45         uint64_t *offset;
  46 } bam_lidx_t;
  47
  48 KSORT_INIT(my_off, pair64_t, pair64_lt);
  49 KHASH_MAP_INIT_INT(my_i, bam_binlist_t);
  50 KHASH_MAP_INIT_STR(s, int)
  51
  52 struct __bam_index_t
  53 {
  54   int32_t n;
  55   khash_t(my_i) **index;
  56   bam_lidx_t *index2;
  57 };
  58
  59 typedef struct __linkbuf_t {
  60         bam1_t b;
  61         uint32_t beg, end;
  62         struct __linkbuf_t *next;
  63 } lbnode_t;
  64
  65 typedef struct {
  66         int cnt, n, max;
  67         lbnode_t **buf;
  68 } mempool_t;
  69
  70 struct __bam_plbuf_t {
  71         mempool_t *mp;
  72         lbnode_t *head, *tail, *dummy;
  73         bam_pileup_f func;
  74         void *func_data;
  75         int32_t tid, pos, max_tid, max_pos;
  76         int max_pu, is_eof;
  77         bam_pileup1_t *pu;
  78         int flag_mask;
  79 };
  80
  81 static mempool_t *mp_init()
  82 {
  83         mempool_t *mp;
  84         mp = (mempool_t*)calloc(1, sizeof(mempool_t));
  85         return mp;
  86 }
  87 static void mp_destroy(mempool_t *mp)
  88 {
  89         int k;
  90         for (k = 0; k < mp->n; ++k) {
  91                 free(mp->buf[k]->b.data);
  92                 free(mp->buf[k]);
  93         }
  94         free(mp->buf);
  95         free(mp);
  96 }
  97 static inline lbnode_t *mp_alloc(mempool_t *mp)
  98 {
  99         ++mp->cnt;
 100         if (mp->n == 0) return (lbnode_t*)calloc(1, sizeof(lbnode_t));
 101         else return mp->buf[--mp->n];
 102 }
 103 static inline void mp_free(mempool_t *mp, lbnode_t *p)
 104 {
 105         --mp->cnt; p->next = 0; // clear lbnode_t::next here
 106         if (mp->n == mp->max) {
 107                 mp->max = mp->max? mp->max<<1 : 256;
 108                 mp->buf = (lbnode_t**)realloc(mp->buf, sizeof(lbnode_t*) * mp->max);
 109         }
 110         mp->buf[mp->n++] = p;
 111 }
 112
 113 static inline int resolve_cigar(bam_pileup1_t *p, uint32_t pos)
 114 {
 115         unsigned k;
 116         bam1_t *b = p->b;
 117         bam1_core_t *c = &b->core;
 118         uint32_t x = c->pos, y = 0;
 119         int ret = 1, is_restart = 1;
 120
 121         if (c->flag&BAM_FUNMAP) return 0; // unmapped read
 122         assert(x <= pos); // otherwise a bug
 123         p->qpos = -1; p->indel = 0; p->is_del = p->is_head = p->is_tail = 0;
 124         for (k = 0; k < c->n_cigar; ++k) {
 125                 int op = bam1_cigar(b)[k] & BAM_CIGAR_MASK; // operation
 126                 int l = bam1_cigar(b)[k] >> BAM_CIGAR_SHIFT; // length
 127                 if (op == BAM_CMATCH) { // NOTE: this assumes the first and the last operation MUST BE a match or a clip
 128                         if (x + l > pos) { // overlap with pos
 129                                 p->indel = p->is_del = 0;
 130                                 p->qpos = y + (pos - x);
 131                                 if (x == pos && is_restart) p->is_head = 1;
 132                                 if (x + l - 1 == pos) { // come to the end of a match
 133                                         if (k < c->n_cigar - 1) { // there are additional operation(s)
 134                                                 uint32_t cigar = bam1_cigar(b)[k+1]; // next CIGAR
 135                                                 int op_next = cigar&BAM_CIGAR_MASK; // next CIGAR operation
 136                                                 if (op_next == BAM_CDEL) p->indel = -(int32_t)(cigar>>BAM_CIGAR_SHIFT); // del
 137                                                 else if (op_next == BAM_CINS) p->indel = cigar>>BAM_CIGAR_SHIFT; // ins
 138                                                 if (op_next == BAM_CSOFT_CLIP || op_next == BAM_CREF_SKIP || op_next == BAM_CHARD_CLIP)
 139                                                         p->is_tail = 1; // tail
 140                                         } else p->is_tail = 1; // this is the last operation; set tail
 141                                 }
 142                         }
 143                         x += l; y += l;
 144                 } else if (op == BAM_CDEL) { // then set ->is_del
 145                         if (x + l > pos) {
 146                                 p->indel = 0; p->is_del = 1;
 147                                 p->qpos = y + (pos - x);
 148                         }
 149                         x += l;
 150                 } else if (op == BAM_CREF_SKIP) x += l;
 151                 else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l;
 152                 is_restart = (op == BAM_CREF_SKIP || op == BAM_CSOFT_CLIP || op == BAM_CHARD_CLIP);
 153                 if (x > pos) {
 154                         if (op == BAM_CREF_SKIP) ret = 0; // then do not put it into pileup at all
 155                         break;
 156                 }
 157         }
 158         assert(x > pos); // otherwise a bug
 159         return ret;
 160
 161 }
 162 // the following code has been taken from bam_plbuf_push
 163 // and modified such that instead of a function call
 164 // the function returns and will continue (if cont is true).
 165 // from where it left off.
 166
 167 // returns
 168 // 1: if buf is full and can be emitted
 169 // 0: if b has been added
 170 // -1: if there was an error
 171 int pysam_pileup_next(const bam1_t *b,
 172                       bam_plbuf_t *buf,
 173                       bam_pileup1_t ** plp,
 174                       int * tid,
 175                       int * pos,
 176                       int * n_plp )
 177 {
 178   *plp = bam_plp_next(buf->iter, tid, pos, n_plp);
 179   if (plp == NULL) return 0;
 180   return 1;
 181 }
 182
 183 typedef struct __bmc_aux_t {
 184         int max;
 185         uint32_t *info;
 186         uint16_t *info16;
 187         errmod_t *em;
 188 } bmc_aux_t;
 189
 190 uint32_t pysam_glf_depth( glf1_t * g )
 191 {
 192   return g->depth;
 193 }
 194
 195
 196 void pysam_dump_glf( glf1_t * g, bam_maqcns_t * c )
 197 {
 198   int x = 0;
 199   fprintf(stderr,
 200           "glf: ref_base=%i, max_mapQ=%i, min_lk=%i, depth=%i",
 201           g->ref_base,
 202           g->max_mapQ,
 203           g->min_lk,
 204           g->depth );
 205
 206   for (x = 0; x < 10; ++x)
 207     fprintf(stderr, ", lk%x=%i, ", x, g->lk[x]);
 208
 209   fprintf(stderr,
 210           "maqcns: het_rate=%f, theta=%f, n_hap=%i, cap_mapQ=%i, errmod=%i, min_baseQ=%i, eta=%f, q_r=%f, aux_max=%i",
 211           c->het_rate,
 212           c->theta,
 213           c->n_hap,
 214           c->cap_mapQ,
 215           c->errmod,
 216           c->min_baseQ,
 217           c->eta,
 218           c->q_r,
 219           c->aux->max);
 220
 221   for (x = 0; x < c->aux->max; ++x)
 222     {
 223       fprintf(stderr, ", info-%i=%i ", x, c->aux->info[x]);
 224       if (c->aux->info[x] == 0) break;
 225     }
 226
 227   for (x = 0; x < c->aux->max; ++x)
 228     {
 229       fprintf(stderr, ", info16-%i=%i ", x, c->aux->info16[x]);
 230       if (c->aux->info16[x] == 0) break;
 231     }
 232 }
 233
 234
 235
 236
 237 // pysam dispatch function to emulate the samtools
 238 // command line within python.
 239 // taken from the main function in bamtk.c
 240 // added code to reset getopt
 241 extern int main_samview(int argc, char *argv[]);
 242 extern int main_import(int argc, char *argv[]);
 243 extern int bam_pileup(int argc, char *argv[]);
 244 extern int bam_merge(int argc, char *argv[]);
 245 extern int bam_sort(int argc, char *argv[]);
 246 extern int bam_index(int argc, char *argv[]);
 247 extern int faidx_main(int argc, char *argv[]);
 248 extern int bam_mating(int argc, char *argv[]);
 249 extern int bam_rmdup(int argc, char *argv[]);
 250 extern int glf3_view_main(int argc, char *argv[]);
 251 extern int bam_flagstat(int argc, char *argv[]);
 252 extern int bam_fillmd(int argc, char *argv[]);
 253
 254 int pysam_dispatch(int argc, char *argv[] )
 255 {
 256   extern int optind;
 257 #ifdef _WIN32
 258   setmode(fileno(stdout), O_BINARY);
 259   setmode(fileno(stdin),  O_BINARY);
 260 #ifdef _USE_KNETFILE
 261   knet_win32_init();
 262 #endif
 263 #endif
 264
 265   // reset getop
 266   optind = 1;
 267
 268   if (argc < 2) return 1;
 269
 270   if (strcmp(argv[1], "view") == 0) return main_samview(argc-1, argv+1);
 271   else if (strcmp(argv[1], "import") == 0) return main_import(argc-1, argv+1);
 272   else if (strcmp(argv[1], "pileup") == 0) return bam_pileup(argc-1, argv+1);
 273   else if (strcmp(argv[1], "merge") == 0) return bam_merge(argc-1, argv+1);
 274   else if (strcmp(argv[1], "sort") == 0) return bam_sort(argc-1, argv+1);
 275   else if (strcmp(argv[1], "index") == 0) return bam_index(argc-1, argv+1);
 276   else if (strcmp(argv[1], "faidx") == 0) return faidx_main(argc-1, argv+1);
 277   else if (strcmp(argv[1], "fixmate") == 0) return bam_mating(argc-1, argv+1);
 278   else if (strcmp(argv[1], "rmdup") == 0) return bam_rmdup(argc-1, argv+1);
 279   else if (strcmp(argv[1], "glfview") == 0) return glf3_view_main(argc-1, argv+1);
 280   else if (strcmp(argv[1], "flagstat") == 0) return bam_flagstat(argc-1, argv+1);
 281   else if (strcmp(argv[1], "calmd") == 0) return bam_fillmd(argc-1, argv+1);
 282   else if (strcmp(argv[1], "fillmd") == 0) return bam_fillmd(argc-1, argv+1);
 283
 284 #if _CURSES_LIB != 0
 285   else if (strcmp(argv[1], "tview") == 0) return bam_tview_main(argc-1, argv+1);
 286 #endif
 287   else
 288     {
 289       fprintf(stderr, "[main] unrecognized command '%s'\n", argv[1]);
 290       return 1;
 291     }
 292   return 0;
 293 }
 294
 295 // taken from samtools/bam_import.c
 296 static inline uint8_t *alloc_data(bam1_t *b, size_t size)
 297 {
 298   if (b->m_data < size)
 299     {
 300       b->m_data = size;
 301       kroundup32(b->m_data);
 302       b->data = (uint8_t*)realloc(b->data, b->m_data);
 303     }
 304   return b->data;
 305 }
 306
 307 // update the variable length data within a bam1_t entry.
 308 // Adds *nbytes_new* - *nbytes_old* into the variable length data of *src* at *pos*.
 309 // Data within the bam1_t entry is moved so that it is
 310 // consistent with the data field lengths.
 311 bam1_t * pysam_bam_update( bam1_t * b,
 312                            const size_t nbytes_old,
 313                            const size_t nbytes_new,
 314                            uint8_t * pos )
 315 {
 316   int d = nbytes_new-nbytes_old;
 317   int new_size;
 318   size_t offset;
 319
 320   // no change
 321   if (d == 0) return b;
 322
 323   new_size = d + b->data_len;
 324   offset = pos - b->data;
 325
 326   //printf("d=%i, old=%i, new=%i, old_size=%i, new_size=%i\n",
 327   // d, nbytes_old, nbytes_new, b->data_len, new_size);
 328
 329   // increase memory if required
 330   if (d > 0)
 331     {
 332       alloc_data( b, new_size );
 333       pos = b->data + offset;
 334     }
 335
 336   if (b->data_len != 0)
 337     {
 338       if (offset < 0 || offset > b->data_len)
 339         fprintf(stderr, "[pysam_bam_insert] illegal offset: '%i'\n", (int)offset);
 340     }
 341
 342   // printf("dest=%p, src=%p, n=%i\n", pos+nbytes_new, pos + nbytes_old, b->data_len - (offset+nbytes_old));
 343   memmove( pos + nbytes_new,
 344            pos + nbytes_old,
 345            b->data_len - (offset + nbytes_old));
 346
 347   b->data_len = new_size;
 348
 349   return b;
 350 }
 351
 352 // translate a nucleotide character to binary code
 353 unsigned char pysam_translate_sequence( const unsigned char s )
 354 {
 355   return bam_nt16_table[s];
 356 }
 357
 358
 359 void bam_init_header_hash(bam_header_t *header);
 360
 361 // translate a reference string *s* to a tid
 362 // code taken from bam_parse_region
 363 int pysam_reference2tid( bam_header_t *header, const char * s )
 364 {
 365
 366   khiter_t iter;
 367   khash_t(s) *h;
 368
 369   bam_init_header_hash(header);
 370   h = (khash_t(s)*)header->hash;
 371
 372   iter = kh_get(s, h, s); /* get the ref_id */
 373   if (iter == kh_end(h)) { // name not found
 374     return -1;
 375   }
 376
 377   return kh_value(h, iter);
 378 }
 379
 380
 381
 382
 383
 384
 385
 386