integrated topmed3 changes from Hyun into main branch

Makefile

@@ -53,6 +53,9 @@ SOURCES = align\
 		index\
 		interval_tree\
 		interval\
+		joint_genotype_sequential\
+		joint_genotyping_buffered_reader\
+		joint_genotyping_record\
 		lfhmm\
 		lhmm\
 		lhmm1\

consolidate_vntrs.cpp

@@ -97,10 +97,13 @@ class Igor : Program
     void consolidate_vntrs()
     {
         bcf1_t *v = vc->odw->get_bcf1_from_pool();
-
+        bcf_hdr_t *h = vc->odr->hdr;
+
         Variant* variant;
         while (vc->odr->read(v))
         {
+            bcf_print_liten(h,v);
+
             variant = new Variant(vc->odw->hdr, v);
             vc->flush_variant_buffer(variant);
             vc->insert_variant_record_into_buffer(variant);

hts_utils.cpp

@@ -112,6 +112,42 @@ void bam_hdr_transfer_contigs_to_bcf_hdr(const bam_hdr_t *sh, bcf_hdr_t *vh)
     if (s.m) free(s.s);
 }
 
+/**
+ * Gets sample names from bam header.
+ */
+std::string bam_hdr_get_sample_name(bam_hdr_t* hdr) {
+  if ( !hdr )
+    error("Failed to read the BAM header");
+
+  const char *p = hdr->text;
+  const char *q, *r;
+  int n = 0;
+  std::string sm;
+  while( ( q = strstr(p, "@RG" ) ) != 0 ) {
+    p = q + 3;
+    r = q = 0;
+    if ( ( q = strstr(p, "\tID:" ) ) != 0 ) q += 4;
+    if ( ( r = strstr(p, "\tSM:" ) ) != 0 ) r += 4;
+    if ( r && q ) {
+      char *u, *v;
+      for (u = (char*)q; *u && *u != '\t' && *u != '\n'; ++u);
+      for (v = (char*)r; *v && *v != '\t' && *v != '\n'; ++v);
+      *u = *v = '\0';
+      if ( sm.empty() )
+	sm = r;
+      else if ( sm.compare(r) != 0 )
+	error("Multiple sample IDs are included in one BAM file - %s, %s", sm.c_str(), r);
+    }
+    else break;
+    p = q > r ? q : r;
+    ++n;
+  }
+  if ( sm.empty() )
+    error("Sample ID information cannot be found");
+  return sm;
+}
+
+
 /**********
  *BAM UTILS
  **********/
@@ -366,10 +402,10 @@ void bam_print(bam_hdr_t *h, bam1_t *s)
 void bcf_hdr_print(bcf_hdr_t *h)
 {
     if ( h->dirty ) bcf_hdr_sync(h);
-    int hlen;
-    char *htxt = bcf_hdr_fmt_text(h, 1, &hlen);
-    std::cerr << htxt;
-    free(htxt);
+    kstring_t htxt = {0,0,0};
+    bcf_hdr_format(h, 0, &htxt);
+    std::cerr << htxt.s;
+    if (htxt.m) free(htxt.s);
 }
 
 /**
@@ -1201,4 +1237,44 @@ float bcf_get_info_flt(bcf_hdr_t *h, bcf1_t *v, const char* tag, float default_v
 int32_t bcf_set_info_flt(bcf_hdr_t *h, bcf1_t *v, const char* tag, float value)
 {
     return bcf_update_info_float(h, v, tag, &value, 1);
+}
+
+/**
+ * Prints a message to STDERR and abort.
+ */
+void error(const char * msg, ...)
+{
+    va_list  ap;    
+    va_start(ap, msg);
+
+    fprintf(stderr, "\nFATAL ERROR - \n");
+    vfprintf(stderr, msg, ap);
+    fprintf(stderr, "\n\n");
+
+    va_end(ap);
+
+    abort();
+    //throw pexception;
+    //exit(EXIT_FAILURE);
+}
+
+/**
+ * Prints a message to STDERR.
+ */
+void notice(const char * msg, ...) 
+{
+    va_list ap;
+    va_start(ap, msg);
+
+    time_t current_time;
+    char buff[255];
+    current_time = time(NULL);
+
+    strftime(buff, 120, "%Y/%m/%d %H:%M:%S", localtime(&current_time));
+
+    fprintf(stderr,"NOTICE [%s] - ", buff);
+    vfprintf(stderr, msg, ap);
+    fprintf(stderr,"\n");
+
+    va_end(ap);
 }

hts_utils.h

@@ -77,6 +77,11 @@ bool str_ends_with(std::string& file_name, const char* ext);
  */
 void bam_hdr_transfer_contigs_to_bcf_hdr(const bam_hdr_t *sh, bcf_hdr_t *vh);
 
+/**
+ * Gets sample names from bam header.
+ */
+std::string bam_hdr_get_sample_name(bam_hdr_t* hdr);
+
 /**
  * Get number of sequences.
  */
@@ -225,7 +230,7 @@ void bcf_hdr_print(bcf_hdr_t *h);
  * Copies contigs found in bcf header to another bcf header.
  */
 void bcf_hdr_transfer_contigs(const bcf_hdr_t *sh, bcf_hdr_t *vh);
-
+    
 /**
  * Checks if a particular header type exists
  * @hdr  - header
@@ -555,4 +560,15 @@ int32_t bcf_set_info_flt(bcf_hdr_t *h, bcf1_t *v, const char* tag, float value);
  */
 #define bcf_set_qual(v, q) ((v)->qual = (q))
 
+
+/**
+ * Prints a message to STDERR and abort.
+ */
+void error(const char * msg, ...);
+
+/**
+ * Gives a notice to STDERR.
+ */
+void notice(const char * msg, ...);
+
 #endif