X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=pysam.git;a=blobdiff_plain;f=pysam%2Fpysam_util.c;fp=pysam%2Fpysam_util.c;h=c9d4543eea72bd3aa193dfff3c7c0f9f5edfac91;hp=91b6fa7e17e04c8fe67a0e68f63292f390f774d4;hb=d02fe5283ed7a93a2f76a5d6dc6e37b40c11b9b1;hpb=d828f9c9aa78e3d1687265b52de841f3f3852089

diff --git a/pysam/pysam_util.c b/pysam/pysam_util.c
index 91b6fa7..c9d4543 100644
--- a/pysam/pysam_util.c
+++ b/pysam/pysam_util.c
@@ -6,6 +6,21 @@
 #include "bam_endian.h"
 #include "knetfile.h"
 #include "pysam_util.h"
+#include "errmod.h" // for pysam_dump 
+
+#ifndef inline
+#define inline __inline
+#endif
+
+// Definition of pysamerr
+#include "stdio.h"
+FILE * pysamerr = NULL;
+
+FILE * pysam_set_stderr( FILE * f )
+{
+  pysamerr = f;
+  return f;
+}
 
 // #######################################################
 // utility routines to avoid using callbacks in bam_fetch
@@ -32,6 +47,7 @@ typedef struct {
 
 KSORT_INIT(my_off, pair64_t, pair64_lt);
 KHASH_MAP_INIT_INT(my_i, bam_binlist_t);
+KHASH_MAP_INIT_STR(s, int)
 
 struct __bam_index_t
 {
@@ -164,6 +180,60 @@ int pysam_pileup_next(const bam1_t *b,
   return 1;
 }
 
+typedef struct __bmc_aux_t {
+	int max;
+	uint32_t *info;
+	uint16_t *info16;
+	errmod_t *em;
+} bmc_aux_t;
+
+uint32_t pysam_glf_depth( glf1_t * g )
+{
+  return g->depth;
+}
+
+
+void pysam_dump_glf( glf1_t * g, bam_maqcns_t * c )
+{
+  int x = 0;
+  fprintf(stderr,
+	  "glf: ref_base=%i, max_mapQ=%i, min_lk=%i, depth=%i",
+	  g->ref_base,
+	  g->max_mapQ,
+	  g->min_lk,
+	  g->depth );
+
+  for (x = 0; x < 10; ++x) 
+    fprintf(stderr, ", lk%x=%i, ", x, g->lk[x]);
+
+  fprintf(stderr,
+	  "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",
+	  c->het_rate,
+	  c->theta,
+	  c->n_hap,
+	  c->cap_mapQ,
+	  c->errmod,
+	  c->min_baseQ,
+	  c->eta,
+	  c->q_r,
+	  c->aux->max);
+  
+  for (x = 0; x < c->aux->max; ++x)
+    {
+      fprintf(stderr, ", info-%i=%i ", x, c->aux->info[x]);
+      if (c->aux->info[x] == 0) break;
+    }
+  
+  for (x = 0; x < c->aux->max; ++x)
+    {
+      fprintf(stderr, ", info16-%i=%i ", x, c->aux->info16[x]);
+      if (c->aux->info16[x] == 0) break;
+    }
+}
+
+
+  
+
 // pysam dispatch function to emulate the samtools
 // command line within python.
 // taken from the main function in bamtk.c
@@ -183,7 +253,7 @@ extern int bam_fillmd(int argc, char *argv[]);
 
 int pysam_dispatch(int argc, char *argv[] )
 {
-
+  extern int optind;
 #ifdef _WIN32
   setmode(fileno(stdout), O_BINARY);
   setmode(fileno(stdin),  O_BINARY);
@@ -191,8 +261,6 @@ int pysam_dispatch(int argc, char *argv[] )
   knet_win32_init();
 #endif
 #endif
-
-  extern int optind;
   
   // reset getop
   optind = 1;
@@ -204,7 +272,7 @@ int pysam_dispatch(int argc, char *argv[] )
   else if (strcmp(argv[1], "pileup") == 0) return bam_pileup(argc-1, argv+1);
   else if (strcmp(argv[1], "merge") == 0) return bam_merge(argc-1, argv+1);
   else if (strcmp(argv[1], "sort") == 0) return bam_sort(argc-1, argv+1);
-  else if (strcmp(argv[1], "index") == 0) return bam_index(argc-1, argv+1);
+  else if (strcmp(argv[1], "index") == 0) return bam_index(argc-1, argv+1); 
   else if (strcmp(argv[1], "faidx") == 0) return faidx_main(argc-1, argv+1);
   else if (strcmp(argv[1], "fixmate") == 0) return bam_mating(argc-1, argv+1);
   else if (strcmp(argv[1], "rmdup") == 0) return bam_rmdup(argc-1, argv+1);
@@ -246,12 +314,14 @@ bam1_t * pysam_bam_update( bam1_t * b,
 			   uint8_t * pos )
 {
   int d = nbytes_new-nbytes_old;
+  int new_size;
+  size_t offset;
 
   // no change
   if (d == 0) return b;
 
-  int new_size = d + b->data_len;
-  size_t offset = pos - b->data;
+  new_size = d + b->data_len;
+  offset = pos - b->data;
 
   //printf("d=%i, old=%i, new=%i, old_size=%i, new_size=%i\n",
   // d, nbytes_old, nbytes_new, b->data_len, new_size);
@@ -286,5 +356,31 @@ unsigned char pysam_translate_sequence( const unsigned char s )
 }
 
 
+void bam_init_header_hash(bam_header_t *header);
+
+// translate a reference string *s* to a tid
+// code taken from bam_parse_region
+int pysam_reference2tid( bam_header_t *header, const char * s )
+{
+  
+  khiter_t iter;
+  khash_t(s) *h;
+  
+  bam_init_header_hash(header);
+  h = (khash_t(s)*)header->hash;
+
+  iter = kh_get(s, h, s); /* get the ref_id */
+  if (iter == kh_end(h)) { // name not found
+    return -1;
+  }
+
+  return kh_value(h, iter);
+}
+
+
+  
+
+
+