[med-svn] [Git][med-team/libseqlib][upstream] New upstream version 1.1.2+git20201210+dfsg

Michael R. Crusoe gitlab at salsa.debian.org
Thu Dec 17 11:46:58 GMT 2020



Michael R. Crusoe pushed to branch upstream at Debian Med / libseqlib


Commits:
438a32ff by Michael R. Crusoe at 2020-12-17T12:16:09+01:00
New upstream version 1.1.2+git20201210+dfsg
- - - - -


28 changed files:

- .gitmodules
- Makefile.am
- Makefile.in
- SeqLib/BWAWrapper.h
- SeqLib/BamReader.h
- SeqLib/BamRecord.h
- SeqLib/BamWalker.h
- SeqLib/BamWriter.h
- SeqLib/FastqReader.h
- SeqLib/FermiAssembler.h
- SeqLib/GenomicRegion.h
- SeqLib/GenomicRegionCollection.cpp
- SeqLib/ReadFilter.h
- + SeqLib/ThreadPool.h
- SeqLib/UnalignedSequence.h
- autogen.sh
- configure
- configure.ac
- src/BWAWrapper.cpp
- src/BamReader.cpp
- src/BamRecord.cpp
- src/BamWriter.cpp
- src/FastqReader.cpp
- src/FermiAssembler.cpp
- src/GenomicRegion.cpp
- src/Makefile.in
- src/ReadFilter.cpp
- src/seqtools/Makefile.in


Changes:

=====================================
.gitmodules
=====================================
@@ -1,11 +1,9 @@
 [submodule "fermi-lite"]
 	path = fermi-lite
-	url = https://github.com/jwalabroad/fermi-lite
+	url = https://github.com/walaj/fermi-lite
 [submodule "htslib"]
 	path = htslib
 	url = https://github.com/samtools/htslib
-        branch = master
 [submodule "bwa"]
 	path = bwa
-	url = https://github.com/jwalabroad/bwa
-	branch = Apache2
+	url = https://github.com/walaj/bwa


=====================================
Makefile.am
=====================================
@@ -2,7 +2,7 @@ AUTOMAKE_OPTIONS = foreign
 SUBDIRS = bwa htslib fermi-lite src
 
 install:  
-	mkdir -p bin && cp src/libseqlib.a fermi-lite/libfml.a bwa/libbwa.a htslib/libhts.a bin 
+	mkdir -p lib && cp src/libseqlib.a fermi-lite/libfml.a bwa/libbwa.a htslib/libhts.a lib
 
 seqtools:
 	mkdir -p bin && cd src/seqtools && make && mv seqtools ../../bin


=====================================
Makefile.in
=====================================
@@ -1,7 +1,7 @@
-# Makefile.in generated by automake 1.15 from Makefile.am.
+# Makefile.in generated by automake 1.16.1 from Makefile.am.
 # @configure_input@
 
-# Copyright (C) 1994-2014 Free Software Foundation, Inc.
+# Copyright (C) 1994-2018 Free Software Foundation, Inc.
 
 # This Makefile.in is free software; the Free Software Foundation
 # gives unlimited permission to copy and/or distribute it,
@@ -132,7 +132,7 @@ am__recursive_targets = \
   $(RECURSIVE_CLEAN_TARGETS) \
   $(am__extra_recursive_targets)
 AM_RECURSIVE_TARGETS = $(am__recursive_targets:-recursive=) TAGS CTAGS \
-	cscope distdir dist dist-all distcheck
+	cscope distdir distdir-am dist dist-all distcheck
 am__tagged_files = $(HEADERS) $(SOURCES) $(TAGS_FILES) \
 	$(LISP)config.h.in
 # Read a list of newline-separated strings from the standard input,
@@ -320,8 +320,8 @@ Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status
 	    echo ' $(SHELL) ./config.status'; \
 	    $(SHELL) ./config.status;; \
 	  *) \
-	    echo ' cd $(top_builddir) && $(SHELL) ./config.status $@ $(am__depfiles_maybe)'; \
-	    cd $(top_builddir) && $(SHELL) ./config.status $@ $(am__depfiles_maybe);; \
+	    echo ' cd $(top_builddir) && $(SHELL) ./config.status $@ $(am__maybe_remake_depfiles)'; \
+	    cd $(top_builddir) && $(SHELL) ./config.status $@ $(am__maybe_remake_depfiles);; \
 	esac;
 
 $(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES)
@@ -454,7 +454,10 @@ distclean-tags:
 	-rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags
 	-rm -f cscope.out cscope.in.out cscope.po.out cscope.files
 
-distdir: $(DISTFILES)
+distdir: $(BUILT_SOURCES)
+	$(MAKE) $(AM_MAKEFLAGS) distdir-am
+
+distdir-am: $(DISTFILES)
 	$(am__remove_distdir)
 	test -d "$(distdir)" || mkdir "$(distdir)"
 	@srcdirstrip=`echo "$(srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \
@@ -519,7 +522,7 @@ distdir: $(DISTFILES)
 	  ! -type d ! -perm -444 -exec $(install_sh) -c -m a+r {} {} \; \
 	|| chmod -R a+r "$(distdir)"
 dist-gzip: distdir
-	tardir=$(distdir) && $(am__tar) | GZIP=$(GZIP_ENV) gzip -c >$(distdir).tar.gz
+	tardir=$(distdir) && $(am__tar) | eval GZIP= gzip $(GZIP_ENV) -c >$(distdir).tar.gz
 	$(am__post_remove_distdir)
 
 dist-bzip2: distdir
@@ -545,7 +548,7 @@ dist-shar: distdir
 	@echo WARNING: "Support for shar distribution archives is" \
 	               "deprecated." >&2
 	@echo WARNING: "It will be removed altogether in Automake 2.0" >&2
-	shar $(distdir) | GZIP=$(GZIP_ENV) gzip -c >$(distdir).shar.gz
+	shar $(distdir) | eval GZIP= gzip $(GZIP_ENV) -c >$(distdir).shar.gz
 	$(am__post_remove_distdir)
 
 dist-zip: distdir
@@ -563,7 +566,7 @@ dist dist-all:
 distcheck: dist
 	case '$(DIST_ARCHIVES)' in \
 	*.tar.gz*) \
-	  GZIP=$(GZIP_ENV) gzip -dc $(distdir).tar.gz | $(am__untar) ;;\
+	  eval GZIP= gzip $(GZIP_ENV) -dc $(distdir).tar.gz | $(am__untar) ;;\
 	*.tar.bz2*) \
 	  bzip2 -dc $(distdir).tar.bz2 | $(am__untar) ;;\
 	*.tar.lz*) \
@@ -573,7 +576,7 @@ distcheck: dist
 	*.tar.Z*) \
 	  uncompress -c $(distdir).tar.Z | $(am__untar) ;;\
 	*.shar.gz*) \
-	  GZIP=$(GZIP_ENV) gzip -dc $(distdir).shar.gz | unshar ;;\
+	  eval GZIP= gzip $(GZIP_ENV) -dc $(distdir).shar.gz | unshar ;;\
 	*.zip*) \
 	  unzip $(distdir).zip ;;\
 	esac
@@ -767,7 +770,7 @@ uninstall-am:
 
 
 install:  
-	mkdir -p bin && cp src/libseqlib.a fermi-lite/libfml.a bwa/libbwa.a htslib/libhts.a bin 
+	mkdir -p lib && cp src/libseqlib.a fermi-lite/libfml.a bwa/libbwa.a htslib/libhts.a lib
 
 seqtools:
 	mkdir -p bin && cd src/seqtools && make && mv seqtools ../../bin


=====================================
SeqLib/BWAWrapper.h
=====================================
@@ -33,6 +33,7 @@ class BWAWrapper {
    */
   BWAWrapper() { 
     idx = 0;
+    copy_comment = false;
     memopt = mem_opt_init();
     memopt->flag |= MEM_F_SOFTCLIP;
   }
@@ -65,6 +66,16 @@ class BWAWrapper {
   void AlignSequence(const std::string& seq, const std::string& name, BamRecordVector& vec, bool hardclip, 
 			   double keep_sec_with_frac_of_primary_score, int max_secondary) const;
 
+  /** Perform a BWA-MEM alignment of a single sequence, and store hits in BamReadVector
+   * @param us UnalignedSequence to be aligned
+   * @param vec Alignment hits are appended to vec
+   * @param hardclip Should the output BamRecord objects be hardclipped
+   * @param keep_sec_with_frac_of_primary_score Set a threshold for whether a secondary alignment should be output
+   * @param max_secondary Set a hard-limit on the number of secondary hits that will be reported
+   */
+  void AlignSequence(const UnalignedSequence& us, BamRecordVector& vec, bool hardclip,
+          double keep_sec_with_frac_of_primary_score, int max_secondary) const;
+
   /** Construct a new bwa index for this object. 
    * @param v vector of references to input (e.g. v = {{"r1", "AT"}};)
    * 
@@ -157,6 +168,11 @@ class BWAWrapper {
 
   /** Check if the index is empty */
   bool IsEmpty() const { return !idx; }
+
+  /** Append FASTA/Q comment to SAM output */
+  void AppendFqComment(){
+      copy_comment = true;
+  }
   
  private:
 
@@ -166,6 +182,9 @@ class BWAWrapper {
   // Store the options in memory
   mem_opt_t * memopt;
 
+  // Store append FASTA/Q comment to SAM output
+  bool copy_comment;
+
   // hold the full index structure
   bwaidx_t* idx;
 


=====================================
SeqLib/BamReader.h
=====================================
@@ -4,16 +4,18 @@
 #include <cassert>
 #include "SeqLib/ReadFilter.h"
 #include "SeqLib/BamWalker.h"
+#include "SeqLib/ThreadPool.h"
 
 // forward declare this from hts.c
 extern "C" {
 int hts_useek(htsFile *file, long uoffset, int where);
 }
 
-class BamReader;
 
 namespace SeqLib {
 
+  class BamReader;
+
   typedef SeqPointer<hts_idx_t> SharedIndex; ///< Shared pointer to the HTSlib index struct
 
   typedef SeqPointer<htsFile> SharedHTSFile; ///< Shared pointer to the HTSlib file pointer
@@ -29,6 +31,11 @@ namespace SeqLib {
 
   _Bam() : m_region_idx(0), empty(true), mark_for_closure(false) {}
 
+    //! Return the header for this BAM
+    const BamHeader& GetHeader() const {
+      return m_hdr;
+    }
+    
     ~_Bam() {}
 
     std::string GetFileName() const { return m_in; }
@@ -45,6 +52,11 @@ namespace SeqLib {
     // the return value here is just passed along from sam_read1
     int32_t load_read(BamRecord& r);
 
+    void set_pool(ThreadPool t) {
+      if (t.IsOpen() && fp) // probably dont need this, it can handle null
+	hts_set_opt(fp.get(),  HTS_OPT_THREAD_POOL, &t.p); //t.p is htsThreadPool
+    }
+
     void reset() {
       empty = true;
       mark_for_closure = false;
@@ -97,7 +109,7 @@ namespace SeqLib {
     bool mark_for_closure;
     
     // open the file pointer
-    bool open_BAM_for_reading();
+    bool open_BAM_for_reading(SeqLib::ThreadPool t);
 
     // hold the reference for CRAM reading
     std::string m_cram_reference;
@@ -141,6 +153,16 @@ class BamReader {
    */
   bool SetRegion(const GenomicRegion& gp);
 
+  /** Assign this BamReader a thread pool
+   * 
+   * The thread pool with stay with this object, but
+   * will not be created or destroyed. This must be done
+   * separately, which allows for multiple readers/writers
+   * to be connected to one thread pool
+   * @return false if the thread pool has not been opened 
+   */
+  bool SetThreadPool(ThreadPool p);
+  
   /** Set up multiple regions. Overwrites current regions. 
    * 
    * This will set the BAM pointer to the first element of the
@@ -271,6 +293,9 @@ class BamReader {
   // hold the reference for CRAM reading
   std::string m_cram_reference;
 
+  // for multicore reading/writing
+  ThreadPool pool;
+
 };
 
 


=====================================
SeqLib/BamRecord.h
=====================================
@@ -286,6 +286,9 @@ class BamRecord {
   /** BamRecord is failed QC */
   inline bool QCFailFlag() const { return b ? ((b->core.flag&BAM_FQCFAIL) != 0) : false; }
 
+  /** BamRecord is supplementary alignment */
+  inline bool SupplementaryFlag() const { return b ? ((b->core.flag&BAM_FSUPPLEMENTARY) != 0) : false; }
+
   /** BamRecord is mapped */
   inline bool MappedFlag() const { return b ? ((b->core.flag&BAM_FUNMAP) == 0) : false; }
 
@@ -332,6 +335,9 @@ class BamRecord {
 
   /** Get the alignment position */
   inline int32_t Position() const { return b ? b->core.pos : -1; }
+
+  /** Get the alignment position, including soft clips */
+  int32_t PositionWithSClips() const;
   
   /** Get the alignment position of mate */
   inline int32_t MatePosition() const { return b ? b->core.mpos: -1; }
@@ -352,6 +358,9 @@ class BamRecord {
   /** Get the end of the alignment */
   int32_t PositionEnd() const;
 
+  /** Get the end of the alignment, including soft clips */
+  int32_t PositionEndWithSClips() const;
+
   /** Get the end of the aligment mate pair */
   int32_t PositionEndMate() const;
 
@@ -364,6 +373,14 @@ class BamRecord {
   /** Get the mapping quality */
   inline int32_t MapQuality() const { return b ? b->core.qual : -1; }
 
+  /** Set the qc fail flag on/off (true -> on) */
+  inline void SetQCFail(bool f) { 
+    if (f)
+      b->core.flag |= BAM_FQCFAIL;
+    else
+      b->core.flag &= ~BAM_FQCFAIL;
+  }
+
   /** Set the mapping quality */
   inline void SetMapQuality(int32_t m) { if (b) b->core.qual = m; }
 
@@ -376,11 +393,21 @@ class BamRecord {
   /** Set the position of the mate read */
   inline void SetPositionMate(int32_t i) { b->core.mpos = i; }
 
-  /** Set the pair mapped flag on */
-  inline void SetPairMappedFlag() { b->core.flag |= BAM_FPAIRED; }
+  /** Set the pair mapped flag on/off (true -> on) */
+  inline void SetPairMappedFlag(bool f) { 
+    if (f)
+      b->core.flag |= BAM_FPAIRED;
+    else
+      b->core.flag &= ~BAM_FPAIRED;
+  }
 
-  /** Set the mate reverse flag on */
-  inline void SetMateReverseFlag() { b->core.flag |= BAM_FMREVERSE; }
+  /** Set the mate reverse flag on/off (true -> on) */
+  inline void SetMateReverseFlag(bool f) { 
+    if (f)
+      b->core.flag |= BAM_FMREVERSE;
+    else
+      b->core.flag &= ~BAM_FMREVERSE;
+  }
 
   /** Get the number of cigar fields */
   inline int32_t CigarSize() const { return b ? b->core.n_cigar : -1; }
@@ -523,7 +550,7 @@ class BamRecord {
   Cigar GetCigar() const {
     uint32_t* c = bam_get_cigar(b);
     Cigar cig;
-    for (int k = 0; k < b->core.n_cigar; ++k) {
+    for (size_t k = 0; k < b->core.n_cigar; ++k) {
       cig.add(CigarField(c[k]));
     }
     return cig;
@@ -596,7 +623,7 @@ class BamRecord {
   inline int32_t AlignmentEndPositionReverse() const {
     uint32_t* c = bam_get_cigar(b);
     int32_t p = 0;
-    for (int32_t i = 0; i < b->core.n_cigar; ++i) { // loop from the end
+    for (size_t i = 0; i < b->core.n_cigar; ++i) { // loop from the end
       if ( (bam_cigar_opchr(c[i]) == 'S') || (bam_cigar_opchr(c[i]) == 'H'))
 	p += bam_cigar_oplen(c[i]);
       else // not a clip, so stop counting
@@ -611,10 +638,10 @@ class BamRecord {
   inline int32_t AlignmentPosition() const {
     uint32_t* c = bam_get_cigar(b);
     int32_t p = 0;
-    for (int32_t i = 0; i < b->core.n_cigar; ++i) {
-      if ( (bam_cigar_opchr(c[i]) == 'S') || (bam_cigar_opchr(c[i]) == 'H'))
+    for (size_t i = 0; i < b->core.n_cigar; ++i) {
+      if (bam_cigar_opchr(c[i]) == 'S')
 	p += bam_cigar_oplen(c[i]);
-      else // not a clip, so stop counting
+      else if (bam_cigar_opchr(c[i]) != 'H') 
 	break;
     }
     return p;
@@ -638,7 +665,7 @@ class BamRecord {
   inline int32_t NumSoftClip() const {
       int32_t p = 0;
       uint32_t* c = bam_get_cigar(b);
-      for (int32_t i = 0; i < b->core.n_cigar; ++i)
+      for (size_t i = 0; i < b->core.n_cigar; ++i)
 	if (bam_cigar_opchr(c[i]) == 'S')
 	  p += bam_cigar_oplen(c[i]);
       return p;
@@ -648,7 +675,7 @@ class BamRecord {
   inline int32_t NumHardClip() const {
       int32_t p = 0;
       uint32_t* c = bam_get_cigar(b);
-      for (int32_t i = 0; i < b->core.n_cigar; ++i) 
+      for (size_t i = 0; i < b->core.n_cigar; ++i) 
 	if (bam_cigar_opchr(c[i]) == 'H')
 	  p += bam_cigar_oplen(c[i]);
       return p;
@@ -659,7 +686,7 @@ class BamRecord {
   inline int32_t NumClip() const {
     int32_t p = 0;
     uint32_t* c = bam_get_cigar(b);
-    for (int32_t i = 0; i < b->core.n_cigar; ++i)
+    for (size_t i = 0; i < b->core.n_cigar; ++i)
       if ( (bam_cigar_opchr(c[i]) == 'S') || (bam_cigar_opchr(c[i]) == 'H') )
 	p += bam_cigar_oplen(c[i]);
     return p;
@@ -715,6 +742,10 @@ class BamRecord {
     if (!p)
       return false;
     t = bam_aux2i(p);
+    int type = *p++;
+    if (!(type == 'i' || type == 'C' || type=='S' || type=='s' || type =='I' || type=='c'))
+      return false;
+
     return true;
   }
 
@@ -727,7 +758,13 @@ class BamRecord {
     uint8_t* p = bam_aux_get(b.get(),tag.c_str());
     if (!p)
       return false;
+
     t = bam_aux2f(p);
+    int type = *p;
+    type = *p++;
+    if (!(type == 'f' || type == 'd'))
+      return false;
+
     return true;
   }
 
@@ -764,7 +801,7 @@ class BamRecord {
   inline std::string CigarString() const {
     std::stringstream cig;
     uint32_t* c = bam_get_cigar(b);
-    for (int k = 0; k < b->core.n_cigar; ++k)
+    for (size_t k = 0; k < b->core.n_cigar; ++k)
       cig << bam_cigar_oplen(c[k]) << "MIDNSHP=XB"[c[k]&BAM_CIGAR_MASK];
     return cig.str();
   }
@@ -792,7 +829,7 @@ class BamRecord {
     if (b->core.tid < 0)
       return std::string();
 
-    if (h.isEmpty())
+    if (!h.isEmpty())
       return h.IDtoName(b->core.tid);
 
     // c++98    


=====================================
SeqLib/BamWalker.h
=====================================
@@ -11,11 +11,6 @@
 #define INT32_MAX 0x7fffffffL
 #endif
 
-extern "C" {
-#include "htslib/cram/cram.h"
-#include "htslib/cram/cram_io.h"
-}
-
 struct idx_delete {
   void operator()(hts_idx_t* x) { if (x) hts_idx_destroy(x); }
 };


=====================================
SeqLib/BamWriter.h
=====================================
@@ -3,6 +3,7 @@
 
 #include <cassert>
 #include "SeqLib/BamRecord.h"
+#include "SeqLib/ThreadPool.h"
 
 namespace SeqLib {
 
@@ -37,6 +38,16 @@ class BamWriter  {
    */
   bool WriteHeader() const;
 
+  /** Assign this BamWriter a thread pool
+   * 
+   * The thread pool with stay with this object, but
+   * will not be created or destroyed. This must be done
+   * separately, which allows for multiple readers/writers
+   * to be connected to one thread pool
+   * @return false if the thread pool has not been opened
+   */
+  bool SetThreadPool(ThreadPool p);
+
   /** Provide a header to this writer 
    * @param h Header for this writer. Copies contents
    */
@@ -93,9 +104,6 @@ class BamWriter  {
   // path to output file
   std::string m_out; 
 
-  // open m_out, true if success
-  void open_BAM_for_writing();
-  
   // output format
   std::string output_format; 
   
@@ -104,6 +112,9 @@ class BamWriter  {
 
   // header
   SeqLib::BamHeader hdr;
+
+  // for multicore reading/writing
+  ThreadPool pool;
   
 };
 


=====================================
SeqLib/FastqReader.h
=====================================
@@ -19,7 +19,10 @@ class FastqReader {
  public:
   
   /** Construct an empty FASTQ/FASTA reader */
- FastqReader() {}
+ FastqReader() {
+   seq = NULL;
+   fp = NULL;
+ }
   
   /** Construct a reader and open a FASTQ/FASTA reader 
    * @param file Path to a FASTQ or FASTA file


=====================================
SeqLib/FermiAssembler.h
=====================================
@@ -13,7 +13,6 @@ extern "C"
 #include "fermi-lite/fml.h"
 #include "fermi-lite/bfc.h"
 }
-
 namespace SeqLib {
 
   /** Sequence assembly using FermiKit from Heng Li
@@ -25,6 +24,8 @@ namespace SeqLib {
     /** Create an empty FermiAssembler with default parameters */
     FermiAssembler ();
 
+    /** Create an Empty FermiAssembler with the provided parameters */
+    FermiAssembler(fml_opt_t &_opt);
     /** Destroy by clearing all stored reads from memory */
     ~FermiAssembler();
 
@@ -83,6 +84,11 @@ namespace SeqLib {
      */
     void SetAggressiveTrim() { opt.mag_opt.flag |= MAG_F_AGGRESSIVE; }
 
+    // Added by Cristian Groza
+    void SetSimplifyBubble() {
+        opt.mag_opt.flag &= ~MAG_F_NO_SIMPL;
+    }
+
     /** From lh3: Drop an overlap if its length is below max_overlap * ratio
      * @param ratio Overlaps below ratio * max_overlap will be removed
      */
@@ -114,6 +120,7 @@ namespace SeqLib {
     /** Return the number of sequences that are controlled by this assembler */
     size_t NumSequences() const { return n_seqs; }
 
+    void WriteGFA(std::ostream &out);
   private:
 
     // reads to assemble


=====================================
SeqLib/GenomicRegion.h
=====================================
@@ -65,9 +65,10 @@ class GenomicRegion {
   GenomicRegion(const std::string& reg, const BamHeader& hdr);
 
   /** Return a string representation of just the first base-pair 
+   * @param h BamHeader used as a lookup table for id to string (i.e. chr name)
    * e.g. 1:10,000
    */
-  std::string PointString() const;
+  std::string PointString(const BamHeader& h) const;
 
   // Randomize the position of this GenomicRegion on the genome
   // 
@@ -94,9 +95,6 @@ class GenomicRegion {
    */
   int32_t DistanceBetweenEnds(const GenomicRegion &gr) const;
 
-  /** Returns identical string as would be obtained from << */
-  std::string ToString() const;
-
   /** Returns true if a.chr < b.chr or a.pos1 < a.pos1 if on same chrome, or if a.pos2 < b.pos2 if same chrom and same pos1 */
   bool operator < (const GenomicRegion& b) const;
 
@@ -129,6 +127,11 @@ class GenomicRegion {
    */
   friend std::ostream& operator<<(std::ostream& out, const GenomicRegion& gr);
 
+  /** Print the chromosome with the correct chromosome name as store in the header 
+   * @param h BamHeader that stores the lookup between chr id and chr name string
+   */
+  std::string ToString(const BamHeader& h) const;
+
   /** Extract the chromosome name as a string 
    * @param h BamHeader to serve as sequence dictionary
    * @exception throws an out_of_range exception if ref id >= h->n_targets


=====================================
SeqLib/GenomicRegionCollection.cpp
=====================================
@@ -9,7 +9,7 @@
 #include <algorithm>
 #include <zlib.h>
 
-#define GZBUFFER 4096
+#define GZBUFFER 65472
 
 //#define DEBUG_OVERLAPS 1
 
@@ -124,7 +124,7 @@ bool GenomicRegionCollection<T>::ReadBED(const std::string & file, const BamHead
 
   gzFile fp = NULL;
   fp = strcmp(file.c_str(), "-")? gzopen(file.c_str(), "r") : gzdopen(fileno(stdin), "r");
-  
+
   if (file.empty() || !fp) {
     std::cerr << "BED file not readable: " << file << std::endl;
     return false;
@@ -215,12 +215,19 @@ bool GenomicRegionCollection<T>::ReadVCF(const std::string & file, const BamHead
     std::string val;
     if (line.empty() || line.at(0) == '#')
       continue;
-    
+
     // read first two columnes
     iss_line >> chr >> pos;
-    
+
     // construct the GenomicRegion
-    T gr(chr, pos, pos, hdr);
+    T gr;
+    try {
+      gr = T(chr, pos, pos, hdr);
+    } catch (...) {
+      std::cerr << "...Could not parse pos: " << pos << std::endl << std::endl
+		<< "...on line " << line << std::endl;
+      
+    }
     if (gr.chr >= 0) 
       m_grv->push_back(gr);
   }


=====================================
SeqLib/ReadFilter.h
=====================================
@@ -549,6 +549,11 @@ class ReadFilterCollection {
     return num;
   }
 
+  /** Check that the filter collection includes something. 
+   * If not, then give it a universal includer
+   */
+  void CheckHasIncluder();
+
   // Return the a tab-delimited tally of which filters were satisfied.
    // Includes the header:
    // total_seen_count total_passed_count region region_passed_count rule rule_passed_count


=====================================
SeqLib/ThreadPool.h
=====================================
@@ -0,0 +1,32 @@
+#ifndef SEQLIB_THREAD_POOL_H
+#define SEQLIB_THREAD_POOL_H
+
+#include <stdexcept>
+#include "SeqLib/BamWalker.h"
+#include "htslib/thread_pool.h"
+
+namespace SeqLib{
+
+class ThreadPool {
+
+ public: 
+  
+ ThreadPool() : nthreads(1) { p.pool = NULL; }
+
+ ThreadPool(int n) : nthreads(1) {
+    p.pool = NULL;
+    if (n < 1)
+      throw std::invalid_argument( "n threads must be > 0");
+    if (!(p.pool = hts_tpool_init(n))) 
+      throw std::runtime_error( "Error creating thread pool");
+  }
+
+  bool IsOpen() { return p.pool != NULL; }
+
+  htsThreadPool p;
+  size_t nthreads;
+
+};
+
+}
+#endif


=====================================
SeqLib/UnalignedSequence.h
=====================================
@@ -16,6 +16,7 @@ extern "C" {
 
 #include <cstring>
 #include <vector>
+#include <iostream>
 
 namespace SeqLib {
 
@@ -30,14 +31,14 @@ namespace SeqLib {
      * @param n Name of the sequence 
      * @param s Sequence, stored as ACTG or N characters
      */
-    UnalignedSequence(const std::string& n, const std::string& s) : Name(n), Seq(s), Qual(std::string()), Strand('*') {}
+    UnalignedSequence(const std::string& n, const std::string& s) : Name(n), Com(std::string()), Seq(s), Qual(std::string()), Strand('*') {}
 
     /** Construct an unaligned sequence with name, sequence and quality score
      * @param n Name of the sequence 
      * @param s Sequence, stored as ACTG or N characters
      * @param q Quality string
      */
-    UnalignedSequence(const std::string& n, const std::string& s, const std::string& q) : Name(n), Seq(s), Qual(q), Strand('*') {}
+    UnalignedSequence(const std::string& n, const std::string& s, const std::string& q) : Name(n), Com(std::string()),  Seq(s), Qual(q), Strand('*') {}
 
     /** Construct an unaligned sequence with name, sequence, quality score and strand
      * @param n Name of the sequence 
@@ -45,12 +46,24 @@ namespace SeqLib {
      * @param q Quality string
      * @param t Strand of the sequence, one of '*', '+', '-'
      */
-    UnalignedSequence(const std::string& n, const std::string& s, const std::string& q, char t) : Name(n), Seq(s), Qual(q), Strand(t) {}
+    UnalignedSequence(const std::string& n, const std::string& s, const std::string& q, char t) : Name(n), Com(std::string()), Seq(s), Qual(q), Strand(t) {}
 
      std::string Name; ///< Name of the contig
+     std::string Com;  ///< Comment of the contig
      std::string Seq;  ///< Sequence of the contig (upper-case ACTGN)
      std::string Qual; ///< Quality scores
      char Strand;      ///< Strand of the sequence. Default is '*'
+     
+    /** Output an unaligned sequence to ostream
+     * @param os ostream
+     * @param us UnalignedSequence
+     */
+    friend std::ostream& operator<<(std::ostream& os, const SeqLib::UnalignedSequence& us){
+        os << "@" << us.Name << " " << us.Com << "\n";
+        os << us.Seq << "\n+\n";
+        os << us.Qual << "\n";
+        return os;
+    }
   };
 
   typedef std::vector<UnalignedSequence> UnalignedSequenceVector; ///< A collection of unaligned sequences


=====================================
autogen.sh
=====================================
@@ -4,4 +4,4 @@ set -ex
 aclocal
 autoconf
 autoheader
-automake -a --add-missing
+automake --add-missing --copy


=====================================
configure
=====================================
@@ -644,7 +644,6 @@ am__nodep
 AMDEPBACKSLASH
 AMDEP_FALSE
 AMDEP_TRUE
-am__quote
 am__include
 DEPDIR
 OBJEXT
@@ -721,7 +720,8 @@ PACKAGE_VERSION
 PACKAGE_TARNAME
 PACKAGE_NAME
 PATH_SEPARATOR
-SHELL'
+SHELL
+am__quote'
 ac_subst_files=''
 ac_user_opts='
 enable_option_checking
@@ -2134,7 +2134,7 @@ ac_link='$CC -o conftest$ac_exeext $CFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $
 ac_compiler_gnu=$ac_cv_c_compiler_gnu
 
 
-am__api_version='1.15'
+am__api_version='1.16'
 
 ac_aux_dir=
 for ac_dir in "$srcdir" "$srcdir/.." "$srcdir/../.."; do
@@ -2679,8 +2679,8 @@ MAKEINFO=${MAKEINFO-"${am_missing_run}makeinfo"}
 
 # For better backward compatibility.  To be removed once Automake 1.9.x
 # dies out for good.  For more background, see:
-# <http://lists.gnu.org/archive/html/automake/2012-07/msg00001.html>
-# <http://lists.gnu.org/archive/html/automake/2012-07/msg00014.html>
+# <https://lists.gnu.org/archive/html/automake/2012-07/msg00001.html>
+# <https://lists.gnu.org/archive/html/automake/2012-07/msg00014.html>
 mkdir_p='$(MKDIR_P)'
 
 # We need awk for the "check" target (and possibly the TAP driver).  The
@@ -2731,7 +2731,7 @@ END
 Aborting the configuration process, to ensure you take notice of the issue.
 
 You can download and install GNU coreutils to get an 'rm' implementation
-that behaves properly: <http://www.gnu.org/software/coreutils/>.
+that behaves properly: <https://www.gnu.org/software/coreutils/>.
 
 If you want to complete the configuration process using your problematic
 'rm' anyway, export the environment variable ACCEPT_INFERIOR_RM_PROGRAM
@@ -3281,45 +3281,45 @@ DEPDIR="${am__leading_dot}deps"
 
 ac_config_commands="$ac_config_commands depfiles"
 
-
-am_make=${MAKE-make}
-cat > confinc << 'END'
+{ $as_echo "$as_me:${as_lineno-$LINENO}: checking whether ${MAKE-make} supports the include directive" >&5
+$as_echo_n "checking whether ${MAKE-make} supports the include directive... " >&6; }
+cat > confinc.mk << 'END'
 am__doit:
-	@echo this is the am__doit target
+	@echo this is the am__doit target >confinc.out
 .PHONY: am__doit
 END
-# If we don't find an include directive, just comment out the code.
-{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for style of include used by $am_make" >&5
-$as_echo_n "checking for style of include used by $am_make... " >&6; }
 am__include="#"
 am__quote=
-_am_result=none
-# First try GNU make style include.
-echo "include confinc" > confmf
-# Ignore all kinds of additional output from 'make'.
-case `$am_make -s -f confmf 2> /dev/null` in #(
-*the\ am__doit\ target*)
-  am__include=include
-  am__quote=
-  _am_result=GNU
-  ;;
-esac
-# Now try BSD make style include.
-if test "$am__include" = "#"; then
-   echo '.include "confinc"' > confmf
-   case `$am_make -s -f confmf 2> /dev/null` in #(
-   *the\ am__doit\ target*)
-     am__include=.include
-     am__quote="\""
-     _am_result=BSD
+# BSD make does it like this.
+echo '.include "confinc.mk" # ignored' > confmf.BSD
+# Other make implementations (GNU, Solaris 10, AIX) do it like this.
+echo 'include confinc.mk # ignored' > confmf.GNU
+_am_result=no
+for s in GNU BSD; do
+  { echo "$as_me:$LINENO: ${MAKE-make} -f confmf.$s && cat confinc.out" >&5
+   (${MAKE-make} -f confmf.$s && cat confinc.out) >&5 2>&5
+   ac_status=$?
+   echo "$as_me:$LINENO: \$? = $ac_status" >&5
+   (exit $ac_status); }
+  case $?:`cat confinc.out 2>/dev/null` in #(
+  '0:this is the am__doit target') :
+    case $s in #(
+  BSD) :
+    am__include='.include' am__quote='"' ;; #(
+  *) :
+    am__include='include' am__quote='' ;;
+esac ;; #(
+  *) :
      ;;
-   esac
-fi
-
-
-{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $_am_result" >&5
-$as_echo "$_am_result" >&6; }
-rm -f confinc confmf
+esac
+  if test "$am__include" != "#"; then
+    _am_result="yes ($s style)"
+    break
+  fi
+done
+rm -f confinc.* confmf.*
+{ $as_echo "$as_me:${as_lineno-$LINENO}: result: ${_am_result}" >&5
+$as_echo "${_am_result}" >&6; }
 
 # Check whether --enable-dependency-tracking was given.
 if test "${enable_dependency_tracking+set}" = set; then :
@@ -4755,6 +4755,122 @@ else
   as_fn_error $? "libz not found, please install zlib (http://www.zlib.net/)" "$LINENO" 5
 fi
 
+{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for library containing lzma_code" >&5
+$as_echo_n "checking for library containing lzma_code... " >&6; }
+if ${ac_cv_search_lzma_code+:} false; then :
+  $as_echo_n "(cached) " >&6
+else
+  ac_func_search_save_LIBS=$LIBS
+cat confdefs.h - <<_ACEOF >conftest.$ac_ext
+/* end confdefs.h.  */
+
+/* Override any GCC internal prototype to avoid an error.
+   Use char because int might match the return type of a GCC
+   builtin and then its argument prototype would still apply.  */
+#ifdef __cplusplus
+extern "C"
+#endif
+char lzma_code ();
+int
+main ()
+{
+return lzma_code ();
+  ;
+  return 0;
+}
+_ACEOF
+for ac_lib in '' lzma; do
+  if test -z "$ac_lib"; then
+    ac_res="none required"
+  else
+    ac_res=-l$ac_lib
+    LIBS="-l$ac_lib  $ac_func_search_save_LIBS"
+  fi
+  if ac_fn_cxx_try_link "$LINENO"; then :
+  ac_cv_search_lzma_code=$ac_res
+fi
+rm -f core conftest.err conftest.$ac_objext \
+    conftest$ac_exeext
+  if ${ac_cv_search_lzma_code+:} false; then :
+  break
+fi
+done
+if ${ac_cv_search_lzma_code+:} false; then :
+
+else
+  ac_cv_search_lzma_code=no
+fi
+rm conftest.$ac_ext
+LIBS=$ac_func_search_save_LIBS
+fi
+{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $ac_cv_search_lzma_code" >&5
+$as_echo "$ac_cv_search_lzma_code" >&6; }
+ac_res=$ac_cv_search_lzma_code
+if test "$ac_res" != no; then :
+  test "$ac_res" = "none required" || LIBS="$ac_res $LIBS"
+
+else
+  as_fn_error $? "liblzma not found, please install lzma (http://www.7-zip.org/sdk.html)" "$LINENO" 5
+fi
+
+{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for library containing BZ2_bzBuffToBuffDecompress" >&5
+$as_echo_n "checking for library containing BZ2_bzBuffToBuffDecompress... " >&6; }
+if ${ac_cv_search_BZ2_bzBuffToBuffDecompress+:} false; then :
+  $as_echo_n "(cached) " >&6
+else
+  ac_func_search_save_LIBS=$LIBS
+cat confdefs.h - <<_ACEOF >conftest.$ac_ext
+/* end confdefs.h.  */
+
+/* Override any GCC internal prototype to avoid an error.
+   Use char because int might match the return type of a GCC
+   builtin and then its argument prototype would still apply.  */
+#ifdef __cplusplus
+extern "C"
+#endif
+char BZ2_bzBuffToBuffDecompress ();
+int
+main ()
+{
+return BZ2_bzBuffToBuffDecompress ();
+  ;
+  return 0;
+}
+_ACEOF
+for ac_lib in '' bz2; do
+  if test -z "$ac_lib"; then
+    ac_res="none required"
+  else
+    ac_res=-l$ac_lib
+    LIBS="-l$ac_lib  $ac_func_search_save_LIBS"
+  fi
+  if ac_fn_cxx_try_link "$LINENO"; then :
+  ac_cv_search_BZ2_bzBuffToBuffDecompress=$ac_res
+fi
+rm -f core conftest.err conftest.$ac_objext \
+    conftest$ac_exeext
+  if ${ac_cv_search_BZ2_bzBuffToBuffDecompress+:} false; then :
+  break
+fi
+done
+if ${ac_cv_search_BZ2_bzBuffToBuffDecompress+:} false; then :
+
+else
+  ac_cv_search_BZ2_bzBuffToBuffDecompress=no
+fi
+rm conftest.$ac_ext
+LIBS=$ac_func_search_save_LIBS
+fi
+{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $ac_cv_search_BZ2_bzBuffToBuffDecompress" >&5
+$as_echo "$ac_cv_search_BZ2_bzBuffToBuffDecompress" >&6; }
+ac_res=$ac_cv_search_BZ2_bzBuffToBuffDecompress
+if test "$ac_res" != no; then :
+  test "$ac_res" = "none required" || LIBS="$ac_res $LIBS"
+
+else
+  as_fn_error $? "lib bz2 not found, please install" "$LINENO" 5
+fi
+
 { $as_echo "$as_me:${as_lineno-$LINENO}: checking for library containing clock_gettime" >&5
 $as_echo_n "checking for library containing clock_gettime... " >&6; }
 if ${ac_cv_search_clock_gettime+:} false; then :
@@ -5567,7 +5683,7 @@ cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
 #
 # INIT-COMMANDS
 #
-AMDEP_TRUE="$AMDEP_TRUE" ac_aux_dir="$ac_aux_dir"
+AMDEP_TRUE="$AMDEP_TRUE" MAKE="${MAKE-make}"
 
 _ACEOF
 
@@ -6181,29 +6297,35 @@ $as_echo "$as_me: executing $ac_file commands" >&6;}
   # Older Autoconf quotes --file arguments for eval, but not when files
   # are listed without --file.  Let's play safe and only enable the eval
   # if we detect the quoting.
-  case $CONFIG_FILES in
-  *\'*) eval set x "$CONFIG_FILES" ;;
-  *)   set x $CONFIG_FILES ;;
-  esac
+  # TODO: see whether this extra hack can be removed once we start
+  # requiring Autoconf 2.70 or later.
+  case $CONFIG_FILES in #(
+  *\'*) :
+    eval set x "$CONFIG_FILES" ;; #(
+  *) :
+    set x $CONFIG_FILES ;; #(
+  *) :
+     ;;
+esac
   shift
-  for mf
+  # Used to flag and report bootstrapping failures.
+  am_rc=0
+  for am_mf
   do
     # Strip MF so we end up with the name of the file.
-    mf=`echo "$mf" | sed -e 's/:.*$//'`
-    # Check whether this is an Automake generated Makefile or not.
-    # We used to match only the files named 'Makefile.in', but
-    # some people rename them; so instead we look at the file content.
-    # Grep'ing the first line is not enough: some people post-process
-    # each Makefile.in and add a new line on top of each file to say so.
-    # Grep'ing the whole file is not good either: AIX grep has a line
+    am_mf=`$as_echo "$am_mf" | sed -e 's/:.*$//'`
+    # Check whether this is an Automake generated Makefile which includes
+    # dependency-tracking related rules and includes.
+    # Grep'ing the whole file directly is not great: AIX grep has a line
     # limit of 2048, but all sed's we know have understand at least 4000.
-    if sed -n 's,^#.*generated by automake.*,X,p' "$mf" | grep X >/dev/null 2>&1; then
-      dirpart=`$as_dirname -- "$mf" ||
-$as_expr X"$mf" : 'X\(.*[^/]\)//*[^/][^/]*/*$' \| \
-	 X"$mf" : 'X\(//\)[^/]' \| \
-	 X"$mf" : 'X\(//\)$' \| \
-	 X"$mf" : 'X\(/\)' \| . 2>/dev/null ||
-$as_echo X"$mf" |
+    sed -n 's,^am--depfiles:.*,X,p' "$am_mf" | grep X >/dev/null 2>&1 \
+      || continue
+    am_dirpart=`$as_dirname -- "$am_mf" ||
+$as_expr X"$am_mf" : 'X\(.*[^/]\)//*[^/][^/]*/*$' \| \
+	 X"$am_mf" : 'X\(//\)[^/]' \| \
+	 X"$am_mf" : 'X\(//\)$' \| \
+	 X"$am_mf" : 'X\(/\)' \| . 2>/dev/null ||
+$as_echo X"$am_mf" |
     sed '/^X\(.*[^/]\)\/\/*[^/][^/]*\/*$/{
 	    s//\1/
 	    q
@@ -6221,53 +6343,48 @@ $as_echo X"$mf" |
 	    q
 	  }
 	  s/.*/./; q'`
-    else
-      continue
-    fi
-    # Extract the definition of DEPDIR, am__include, and am__quote
-    # from the Makefile without running 'make'.
-    DEPDIR=`sed -n 's/^DEPDIR = //p' < "$mf"`
-    test -z "$DEPDIR" && continue
-    am__include=`sed -n 's/^am__include = //p' < "$mf"`
-    test -z "$am__include" && continue
-    am__quote=`sed -n 's/^am__quote = //p' < "$mf"`
-    # Find all dependency output files, they are included files with
-    # $(DEPDIR) in their names.  We invoke sed twice because it is the
-    # simplest approach to changing $(DEPDIR) to its actual value in the
-    # expansion.
-    for file in `sed -n "
-      s/^$am__include $am__quote\(.*(DEPDIR).*\)$am__quote"'$/\1/p' <"$mf" | \
-	 sed -e 's/\$(DEPDIR)/'"$DEPDIR"'/g'`; do
-      # Make sure the directory exists.
-      test -f "$dirpart/$file" && continue
-      fdir=`$as_dirname -- "$file" ||
-$as_expr X"$file" : 'X\(.*[^/]\)//*[^/][^/]*/*$' \| \
-	 X"$file" : 'X\(//\)[^/]' \| \
-	 X"$file" : 'X\(//\)$' \| \
-	 X"$file" : 'X\(/\)' \| . 2>/dev/null ||
-$as_echo X"$file" |
-    sed '/^X\(.*[^/]\)\/\/*[^/][^/]*\/*$/{
-	    s//\1/
-	    q
-	  }
-	  /^X\(\/\/\)[^/].*/{
+    am_filepart=`$as_basename -- "$am_mf" ||
+$as_expr X/"$am_mf" : '.*/\([^/][^/]*\)/*$' \| \
+	 X"$am_mf" : 'X\(//\)$' \| \
+	 X"$am_mf" : 'X\(/\)' \| . 2>/dev/null ||
+$as_echo X/"$am_mf" |
+    sed '/^.*\/\([^/][^/]*\)\/*$/{
 	    s//\1/
 	    q
 	  }
-	  /^X\(\/\/\)$/{
+	  /^X\/\(\/\/\)$/{
 	    s//\1/
 	    q
 	  }
-	  /^X\(\/\).*/{
+	  /^X\/\(\/\).*/{
 	    s//\1/
 	    q
 	  }
 	  s/.*/./; q'`
-      as_dir=$dirpart/$fdir; as_fn_mkdir_p
-      # echo "creating $dirpart/$file"
-      echo '# dummy' > "$dirpart/$file"
-    done
+    { echo "$as_me:$LINENO: cd "$am_dirpart" \
+      && sed -e '/# am--include-marker/d' "$am_filepart" \
+        | $MAKE -f - am--depfiles" >&5
+   (cd "$am_dirpart" \
+      && sed -e '/# am--include-marker/d' "$am_filepart" \
+        | $MAKE -f - am--depfiles) >&5 2>&5
+   ac_status=$?
+   echo "$as_me:$LINENO: \$? = $ac_status" >&5
+   (exit $ac_status); } || am_rc=$?
   done
+  if test $am_rc -ne 0; then
+    { { $as_echo "$as_me:${as_lineno-$LINENO}: error: in \`$ac_pwd':" >&5
+$as_echo "$as_me: error: in \`$ac_pwd':" >&2;}
+as_fn_error $? "Something went wrong bootstrapping makefile fragments
+    for automatic dependency tracking.  Try re-running configure with the
+    '--disable-dependency-tracking' option to at least be able to build
+    the package (albeit without support for automatic dependency tracking).
+See \`config.log' for more details" "$LINENO" 5; }
+  fi
+  { am_dirpart=; unset am_dirpart;}
+  { am_filepart=; unset am_filepart;}
+  { am_mf=; unset am_mf;}
+  { am_rc=; unset am_rc;}
+  rm -f conftest-deps.mk
 }
  ;;
 


=====================================
configure.ac
=====================================
@@ -18,6 +18,8 @@ AC_CHECK_HEADER([zlib.h])
 
 # Check for libraries
 AC_SEARCH_LIBS([gzopen],[z],,[AC_MSG_ERROR([libz not found, please install zlib (http://www.zlib.net/)])])
+AC_SEARCH_LIBS([lzma_code],[lzma],,[AC_MSG_ERROR([liblzma not found, please install lzma (http://www.7-zip.org/sdk.html)])])
+AC_SEARCH_LIBS([BZ2_bzBuffToBuffDecompress],[bz2],,[AC_MSG_ERROR([lib bz2 not found, please install])])
 AC_SEARCH_LIBS([clock_gettime], [rt], [AC_DEFINE([HAVE_CLOCK_GETTIME], [1], [clock_getttime found])], )
 
 # check for c++11


=====================================
src/BWAWrapper.cpp
=====================================
@@ -16,6 +16,25 @@ extern "C" {
   #include <string.h>
 }
 
+// custom sort 
+bool aln_sort(const mem_aln_t& lhs, const mem_aln_t& rhs)
+{
+  if (lhs.mapq > rhs.mapq) //want descending order of mapq
+    return true;
+  else if (lhs.mapq < rhs.mapq)
+    return false;
+  // break tied mapq with position
+  if (lhs.rid < rhs.rid)
+    return true;
+  else if (lhs.rid > rhs.rid)
+    return false;
+  if (lhs.pos < rhs.pos)
+    return true;
+  else if (lhs.pos > rhs.pos)
+    return false;
+  return false;
+}
+
 //#define DEBUG_BWATOOLS 1
 
 #define _set_pac(pac, l, c) ((pac)[(l)>>2] |= (c)<<((~(l)&3)<<1))
@@ -244,13 +263,25 @@ namespace SeqLib {
     ar = mem_align1(memopt, idx->bwt, idx->bns, idx->pac, seq.length(), seq.data()); // get all the hits (was c_str())
 
 #ifdef DEBUG_BWATOOLS
-    std::cout << "num hits: " << ar.n << std::endl;
+      std::cerr << "num hits: " << ar.n << " " << name << std::endl;
+      size_t secondary_count_debug = 0;
+      for (size_t i = 0; i < ar.n; ++i) {
+	if (ar.a[i].secondary >= 0)
+	  secondary_count_debug++;
+      }
+      std::cerr << "secondary count " << secondary_count_debug << std::endl;
+    }
 #endif    
 
     double primary_score = 0;
 
     int secondary_count = 0;
 
+    //setup a vector to store the hits. Why remake a vector? Because
+    //I want them to be in the same order for all runs. With BWA multithreading,
+    //it may not be so in the initial ar.a vector.
+    std::vector<mem_aln_t> a;
+
     //size_t num_secondary = 0;
     // loop through the hits
     for (size_t i = 0; i < ar.n; ++i) {
@@ -259,34 +290,40 @@ namespace SeqLib {
       	continue; // skip secondary alignments
       
       // get forward-strand position and CIGAR
-      mem_aln_t a;
+      mem_aln_t a_aln;
 
-      a = mem_reg2aln(memopt, idx->bns, idx->pac, seq.length(), seq.c_str(), &ar.a[i]); 
+      a_aln = mem_reg2aln(memopt, idx->bns, idx->pac, seq.length(), seq.c_str(), &ar.a[i]); 
 
+      a.push_back(a_aln);
+    }
+    
+    // sort it 
+    std::sort(a.begin(), a.end(), aln_sort);
+    
+    for (size_t i = 0; i < a.size(); ++i) {
+      
       // if score not sufficient or past cap, continue
-      bool sec_and_low_score =  ar.a[i].secondary >= 0 && (primary_score * keep_sec_with_frac_of_primary_score) > a.score;
-      bool sec_and_cap_hit = ar.a[i].secondary >= 0 && (int)i > max_secondary;
+      //bool sec_and_low_score =  ar.a[i].secondary >= 0 && (primary_score * keep_sec_with_frac_of_primary_score) > a.score;
+      //bool sec_and_cap_hit = ar.a[i].secondary >= 0 && (int)i > max_secondary;
+      bool sec_and_low_score = (a[i].flag&BAM_FSECONDARY) && (primary_score * keep_sec_with_frac_of_primary_score) > a[i].score;
+      bool sec_and_cap_hit =   (a[i].flag&BAM_FSECONDARY) && (int)i > max_secondary;
       if (sec_and_low_score || sec_and_cap_hit) {
-	free(a.cigar);
+	free(a[i].cigar);
 	continue;
-      } else if (ar.a[i].secondary < 0) {
-	primary_score = a.score;
+      } else if (!(a[i].flag&BAM_FSECONDARY)) {
+	primary_score = a[i].score;
 	//num_secondary = 0;
       }
 
-#ifdef DEBUG_BWATOOLS
-      std::cerr << "allocing bamread" << std::endl;
-#endif
-
       // instantiate the read
       BamRecord b;
       b.init();
 
-      b.b->core.tid = a.rid;
-      b.b->core.pos = a.pos;
-      b.b->core.qual = a.mapq;
-      b.b->core.flag = a.flag;
-      b.b->core.n_cigar = a.n_cigar;
+      b.b->core.tid = a[i].rid;
+      b.b->core.pos = a[i].pos;
+      b.b->core.qual = a[i].mapq;
+      b.b->core.flag = a[i].flag;
+      b.b->core.n_cigar = a[i].n_cigar;
       
       // set dumy mate
       b.b->core.mtid = -1;
@@ -294,7 +331,7 @@ namespace SeqLib {
       b.b->core.isize = 0;
 
       // if alignment is reverse, set it
-      if (a.is_rev) 
+      if (a[i].is_rev) 
 	b.b->core.flag |= BAM_FREVERSE;
 
       std::string new_seq = seq;
@@ -302,11 +339,11 @@ namespace SeqLib {
       if (hardclip) {
 	size_t tstart = 0;
 	size_t len = 0;
-	for (int i = 0; i < a.n_cigar; ++i) {
-	  if (i == 0 && bam_cigar_op(a.cigar[i]) == BAM_CREF_SKIP) // first N (e.g. 20N50M)
-	    tstart = bam_cigar_oplen(a.cigar[i]);
-	  else if (bam_cigar_type(bam_cigar_op(a.cigar[i]))&1) // consumes query, but not N
-	    len += bam_cigar_oplen(a.cigar[i]);
+	for (int i = 0; i < a[i].n_cigar; ++i) {
+	  if (i == 0 && bam_cigar_op(a[i].cigar[i]) == BAM_CREF_SKIP) // first N (e.g. 20N50M)
+	    tstart = bam_cigar_oplen(a[i].cigar[i]);
+	  else if (bam_cigar_type(bam_cigar_op(a[i].cigar[i]))&1) // consumes query, but not N
+	    len += bam_cigar_oplen(a[i].cigar[i]);
 	}
 	assert(len > 0);
 	assert(tstart + len <= seq.length());
@@ -316,19 +353,15 @@ namespace SeqLib {
       // allocate all the data
       b.b->core.l_qname = name.length() + 1;
       b.b->core.l_qseq = new_seq.length(); //(seq.length()>>1) + seq.length() % 2; // 4-bit encoding
-      b.b->l_data = b.b->core.l_qname + (a.n_cigar<<2) + ((b.b->core.l_qseq+1)>>1) + (b.b->core.l_qseq);
+      b.b->l_data = b.b->core.l_qname + (a[i].n_cigar<<2) + ((b.b->core.l_qseq+1)>>1) + (b.b->core.l_qseq);
       b.b.get()->data = (uint8_t*)malloc(b.b.get()->l_data);
 
-#ifdef DEBUG_BWATOOLS
-      std::cerr << "memcpy" << std::endl;
-#endif
-
       // allocate the qname
       memcpy(b.b->data, name.c_str(), name.length() + 1);
 
       // allocate the cigar. Reverse if aligned to neg strand, since mem_aln_t stores
       // cigars relative to referemce string oreiatnion, not forward alignment
-      memcpy(b.b->data + b.b->core.l_qname, (uint8_t*)a.cigar, a.n_cigar<<2);
+      memcpy(b.b->data + b.b->core.l_qname, (uint8_t*)a[i].cigar, a[i].n_cigar<<2);
 
       // convert N to S or H
       int new_val = hardclip ? BAM_CHARD_CLIP : BAM_CSOFT_CLIP;
@@ -346,7 +379,7 @@ namespace SeqLib {
       // TODO move this out of bigger loop
       int slen = new_seq.length();
       int j = 0;
-      if (a.is_rev) {
+      if (a[i].is_rev) {
 	for (int i = slen-1; i >= 0; --i) {
 	  
 	  // bad idea but works for now
@@ -384,23 +417,19 @@ namespace SeqLib {
 	}
       }
 
-#ifdef DEBUG_BWATOOLS
-      std::cerr << "memcpy3" << std::endl;
-#endif
-
       // allocate the quality to NULL
       uint8_t* s = bam_get_qual(b.b);
       s[0] = 0xff;
 
       b.AddIntTag("NA", ar.n); // number of matches
-      b.AddIntTag("NM", a.NM);
+      b.AddIntTag("NM", a[i].NM);
 
-      if (a.XA)
-	b.AddZTag("XA", std::string(a.XA));
+      if (a[i].XA)
+	b.AddZTag("XA", std::string(a[i].XA));
 
       // add num sub opt
-      b.AddIntTag("SB", ar.a[i].sub_n);
-      b.AddIntTag("AS", a.score);
+      //b.AddIntTag("SB", ar.a[i].sub_n);
+      b.AddIntTag("AS", a[i].score);
 
       // count num secondaries
       if (b.SecondaryFlag())
@@ -410,24 +439,34 @@ namespace SeqLib {
 
 #ifdef DEBUG_BWATOOLS
       // print alignment
-      printf("\t%c\t%s\t%ld\t%d\t", "+-"[a.is_rev], idx->bns->anns[a.rid].name, (long)a.pos, a.mapq);
-      for (int k = 0; k < a.n_cigar; ++k) // print CIGAR
-      	printf("%d%c", a.cigar[k]>>4, "MIDSH"[a.cigar[k]&0xf]);
-      printf("\t%d\n", a.NM); // print edit distance
-      std::cerr << "final done" << std::endl;
+      printf("\t%c\t%s\t%ld\t%d\t", "+-"[a[i].is_rev], idx->bns->anns[a[i].rid].name, (long)a[i].pos, a[i].mapq);
+      for (int k = 0; k < a[i].n_cigar; ++k) // print CIGAR
+      	printf("%d%c", a[i].cigar[k]>>4, "MIDSH"[a[i].cigar[k]&0xf]);
+      printf("\t%d\n", a[i].NM); // print edit distance
 #endif
       
-      free(a.cigar); // don't forget to deallocate CIGAR
+      free(a[i].cigar); // don't forget to deallocate CIGAR
     }
+    
     free (ar.a); // dealloc the hit list
-
+    
     // add the secondary counts
     for (BamRecordVector::iterator i = vec.begin(); i != vec.end(); ++i)
       i->AddIntTag("SQ", secondary_count);
     
+  }
+  
+  void BWAWrapper::AlignSequence(const UnalignedSequence& us, BamRecordVector& vec, bool hardclip,
+          double keep_sec_with_frac_of_primary_score, int max_secondary) const{
+      // call overloaded AlignSequence to do alignment
+      AlignSequence(us.Seq, us.Name, vec, hardclip, keep_sec_with_frac_of_primary_score, max_secondary);
+      // append FASTA/Q comment to SAM output if needed
+      for (BamRecordVector::iterator i = vec.begin(); i != vec.end(); ++i)
+          if (copy_comment)
+              i->AddZTag("BC", us.Com);
 }
-
-  // modified from bwa (heng li)
+  
+// modified from bwa (heng li)
 uint8_t* BWAWrapper::seqlib_add1(const kseq_t *seq, bntseq_t *bns, uint8_t *pac, int64_t *m_pac, int *m_seqs, int *m_holes, bntamb1_t **q)
 {
   bntann1_t *p;


=====================================
src/BamReader.cpp
=====================================
@@ -1,6 +1,5 @@
 #include "SeqLib/BamReader.h"
 
-
 //#define DEBUG_WALKER 1
 
 namespace SeqLib {
@@ -149,8 +148,9 @@ void BamReader::Reset() {
       return false;
     
     _Bam new_bam(bam);
+    if (!m_cram_reference.empty()) new_bam.m_cram_reference = m_cram_reference;
     new_bam.m_region = &m_region;
-    bool success = new_bam.open_BAM_for_reading();
+    bool success = new_bam.open_BAM_for_reading(pool);
     m_bams.insert(std::pair<std::string, _Bam>(bam, new_bam));
     return success;
   }
@@ -173,22 +173,32 @@ BamReader::BamReader() {}
 
   }
 
+  bool BamReader::SetThreadPool(ThreadPool p) {
+    if (!p.IsOpen())
+      return false;
+    pool = p;
+    for (_BamMap::iterator b = m_bams.begin(); b != m_bams.end(); ++b)
+      b->second.set_pool(p);
+    return true;
+  }
+  
   BamHeader BamReader::Header() const { 
     if (m_bams.size()) 
       return m_bams.begin()->second.m_hdr; 
     return BamHeader(); 
   }
 
-  bool _Bam::open_BAM_for_reading() {
+  bool _Bam::open_BAM_for_reading(SeqLib::ThreadPool t) {
 
     // HTS open the reader
     fp = SharedHTSFile(hts_open(m_in.c_str(), "r"), htsFile_delete()); 
 
+    // connect the thread pool (may already be done, but its ok
+    set_pool(t);
+
     // open cram reference
     if (!m_cram_reference.empty()) {
-      char * m_cram_reference_cstr = strdup(m_cram_reference.c_str());
-      int ret = cram_load_reference(fp->fp.cram, m_cram_reference_cstr);
-      free(m_cram_reference_cstr);
+      int ret = hts_set_fai_filename(fp.get(), m_cram_reference.c_str());
       if (ret < 0) 
 	throw std::invalid_argument("Could not read reference genome " + m_cram_reference + " for CRAM opt");
     }


=====================================
src/BamRecord.cpp
=====================================
@@ -41,10 +41,29 @@ namespace SeqLib {
     b = SeqPointer<bam1_t>(a, free_delete()); 
   }
 
+  int32_t BamRecord::PositionWithSClips() const {
+    if(!b) return -1; // to be consistent with BamRecord::Position()
+
+    uint32_t* cig = bam_get_cigar(b);
+    return ((*cig) & 0xF) == BAM_CSOFT_CLIP ? b->core.pos - ((*cig) >> 4) : b->core.pos;
+  }
+
   int32_t BamRecord::PositionEnd() const { 
     return b ? (b->core.l_qseq > 0 ? bam_endpos(b.get()) : b->core.pos + GetCigar().NumQueryConsumed()) : -1;
   }
 
+  int32_t BamRecord::PositionEndWithSClips() const {
+    if(!b) return -1; // to be consistent with BamRecord::PositionEnd()
+
+    uint32_t* cig_last = bam_get_cigar(b) + b->core.n_cigar - 1;
+    if(b->core.l_qseq > 0) {
+      return ((*cig_last) & 0xF) == BAM_CSOFT_CLIP ? bam_endpos(b.get()) + ((*cig_last) >> 4) :
+                                                     bam_endpos(b.get());
+    } else {
+      return b->core.pos + GetCigar().NumQueryConsumed();
+    }
+  }
+
   int32_t BamRecord::PositionEndMate() const { 
     return b ? (b->core.mpos + (b->core.l_qseq > 0 ? b->core.l_qseq : GetCigar().NumQueryConsumed())) : -1;
   }
@@ -453,8 +472,10 @@ namespace SeqLib {
   }
 
   bool BamRecord::GetTag(const std::string& tag, std::string& s) const {
+
     if (GetZTag(tag, s))
       return true;
+
     int32_t t;
     if (GetIntTag(tag, t)) {
       std::stringstream ss; 
@@ -462,6 +483,7 @@ namespace SeqLib {
       s = ss.str();
       return true;
     } 
+
     float f;
     if (GetFloatTag(tag, f)) {
       std::stringstream ss; 
@@ -469,13 +491,20 @@ namespace SeqLib {
       s = ss.str();
       return true;
     }     
+
     return false;
+
   }
 
   bool BamRecord::GetZTag(const std::string& tag, std::string& s) const {
     uint8_t* p = bam_aux_get(b.get(),tag.c_str());
     if (!p)
       return false;
+
+    int type = *p; 
+    if (type != 'Z')
+      return false;
+    
     char* pp = bam_aux2Z(p);
     if (!pp) 
       return false;


=====================================
src/BamWriter.cpp
=====================================
@@ -77,6 +77,9 @@ bool BamWriter::BuildIndex() const {
     // hts open the writer
     fop = SeqPointer<htsFile>(hts_open(m_out.c_str(), output_format.c_str()), htsFile_delete());
 
+    // open the thread pool. It's OK if already connected before opening
+    SetThreadPool(pool);
+
     if (!fop) {
       return false;
       //throw std::runtime_error("BamWriter::Open - Cannot open output file: " + f);
@@ -117,6 +120,14 @@ std::ostream& operator<<(std::ostream& out, const BamWriter& b)
   return out;
 }
 
+bool BamWriter::SetThreadPool(ThreadPool p) {
+  if (!p.IsOpen()) 
+    return false;
+  pool = p;
+  if (fop.get())
+    hts_set_opt(fop.get(),  HTS_OPT_THREAD_POOL, &pool.p);
+  return true;
+}
   //does not return false if file not found
 bool BamWriter::SetCramReference(const std::string& ref) {
 


=====================================
src/FastqReader.cpp
=====================================
@@ -47,6 +47,8 @@ bool FastqReader::GetNextSequence(UnalignedSequence& s) {
 
   if (seq->name.s)
     s.Name = std::string(seq->name.s, seq->name.l);
+  if (seq->comment.s)
+      s.Com = std::string(seq->comment.s, seq->comment.l);
   if (seq->seq.s)
     s.Seq = std::string(seq->seq.s, seq->seq.l);
   if (seq->qual.s)


=====================================
src/FermiAssembler.cpp
=====================================
@@ -1,4 +1,5 @@
 #include "SeqLib/FermiAssembler.h"
+#include "fermi-lite/fml.h"
 #define MAG_MIN_NSR_COEF .1
 
 namespace SeqLib {
@@ -6,171 +7,199 @@ namespace SeqLib {
   FermiAssembler::FermiAssembler()  : m_seqs(0), m(0), size(0), n_seqs(0), n_utg(0), m_utgs(0)  {
     fml_opt_init(&opt);
   }
-  
-  FermiAssembler::~FermiAssembler() {
-    ClearReads();
-    ClearContigs();
-  }
-
-  // code copied and slightly modified from 
-  // fermi-lite/misc.c by Heng Li
-  void FermiAssembler::DirectAssemble(float kcov) {
 
-    rld_t *e = fml_seq2fmi(&opt, n_seqs, m_seqs);
-    mag_t *g = fml_fmi2mag(&opt, e);
-
-    opt.mag_opt.min_ensr = opt.mag_opt.min_ensr > kcov * MAG_MIN_NSR_COEF? opt.mag_opt.min_ensr : (int)(kcov * MAG_MIN_NSR_COEF + .499);
-    //opt.mag_opt.min_ensr = opt.mag_opt.min_ensr < opt0->max_cnt? opt.mag_opt.min_ensr : opt0->max_cnt;
-    //opt.mag_opt.min_ensr = opt.mag_opt.min_ensr > opt0->min_cnt? opt.mag_opt.min_ensr : opt0->min_cnt;
-    opt.mag_opt.min_insr = opt.mag_opt.min_ensr - 1;
-    fml_mag_clean(&opt, g);
-    m_utgs = fml_mag2utg(g, &n_utg);
+  FermiAssembler::FermiAssembler(fml_opt_t &_opt) :
+      m_seqs(0), m(0), size(0), n_seqs(0), n_utg(0), m_utgs(0), opt(_opt)
+  {
   }
 
-  void FermiAssembler::AddRead(const BamRecord& r) {
-    AddRead(UnalignedSequence(r.Qname(), r.Sequence(), r.Qualities())); // probably faster way
-  }
+    FermiAssembler::~FermiAssembler()
+    {
+      ClearReads();
+      ClearContigs();
+    }
 
-  void FermiAssembler::AddRead(const UnalignedSequence& r) {
-
-    if (r.Seq.empty())
-      return;
-    if (r.Name.empty())
-      return;
-
-    // dynamically alloc the memory
-    if (m <= n_seqs)
-      m = m <= 0 ? 32 : (m*2); // if out of mem, double it
-    m_seqs = (fseq1_t*)realloc(m_seqs, m * sizeof(fseq1_t));
-    
-    // add the name
-    m_names.push_back(r.Name);
-
-    // construct the seq
-    fseq1_t *s;
-    s = &m_seqs[n_seqs];
-    s->seq   = strdup(r.Seq.c_str());
-    s->qual = r.Qual.empty() ? NULL : strdup(r.Qual.c_str());
-    
-    s->l_seq = r.Seq.length();
-    size += m_seqs[n_seqs++].l_seq;
+    // code copied and slightly modified from
+    // fermi-lite/misc.c by Heng Li
+    void FermiAssembler::DirectAssemble(float kcov)
+    {
+      rld_t *e = fml_seq2fmi(&opt, n_seqs, m_seqs);
+      mag_t *g = fml_fmi2mag(&opt, e);
+
+      opt.mag_opt.min_ensr = opt.mag_opt.min_ensr > kcov * MAG_MIN_NSR_COEF
+                               ? opt.mag_opt.min_ensr
+                               : (int)(kcov * MAG_MIN_NSR_COEF + .499);
+      // opt.mag_opt.min_ensr = opt.mag_opt.min_ensr < opt0->max_cnt?
+      // opt.mag_opt.min_ensr : opt0->max_cnt; opt.mag_opt.min_ensr =
+      // opt.mag_opt.min_ensr > opt0->min_cnt? opt.mag_opt.min_ensr :
+      // opt0->min_cnt;
+      opt.mag_opt.min_insr = opt.mag_opt.min_ensr - 1;
+      fml_mag_clean(&opt, g);
+      m_utgs = fml_mag2utg(g, &n_utg);
+    }
 
-  }
-  
-  void FermiAssembler::AddReads(const UnalignedSequenceVector& v) {
+    void FermiAssembler::AddRead(const BamRecord &r)
+    {
+      AddRead(UnalignedSequence(
+        r.Qname(), r.Sequence(), r.Qualities())); // probably faster way
+    }
 
-    // alloc the memory
-    m = n_seqs + v.size();
-    m_seqs = (fseq1_t*)realloc(m_seqs, m * sizeof(fseq1_t));
+    void FermiAssembler::AddRead(const UnalignedSequence &r)
+    {
+      if (r.Seq.empty()) return;
+      if (r.Name.empty()) return;
 
-    for (UnalignedSequenceVector::const_iterator r = v.begin(); r != v.end(); ++r) {
-      m_names.push_back(r->Name);
-      fseq1_t *s;
+      // dynamically alloc the memory
+      if (m <= n_seqs) m = m <= 0 ? 32 : (m * 2); // if out of mem, double it
+      m_seqs = (fseq1_t *)realloc(m_seqs, m * sizeof(fseq1_t));
 
-      s = &m_seqs[n_seqs];
+      // add the name
+      m_names.push_back(r.Name);
 
-      s->seq   = strdup(r->Seq.c_str());
-      s->qual  = strdup(r->Qual.c_str());
+      // construct the seq
+      fseq1_t *s;
+      s       = &m_seqs[n_seqs];
+      s->seq  = strdup(r.Seq.c_str());
+      s->qual = r.Qual.empty() ? NULL : strdup(r.Qual.c_str());
 
-      s->l_seq = r->Seq.length();
+      s->l_seq = r.Seq.length();
       size += m_seqs[n_seqs++].l_seq;
     }
 
+    void FermiAssembler::AddReads(const UnalignedSequenceVector &v)
+    {
+      // alloc the memory
+      m      = n_seqs + v.size();
+      m_seqs = (fseq1_t *)realloc(m_seqs, m * sizeof(fseq1_t));
 
-  }
-  void FermiAssembler::AddReads(const BamRecordVector& brv) {
+      for (UnalignedSequenceVector::const_iterator r = v.begin(); r != v.end();
+           ++r) {
+        m_names.push_back(r->Name);
+        fseq1_t *s;
 
-    // alloc the memory
-    m_seqs = (fseq1_t*)realloc(m_seqs, (n_seqs + brv.size()) * sizeof(fseq1_t));
+        s = &m_seqs[n_seqs];
 
-    uint64_t size = 0;
-    for (BamRecordVector::const_iterator r = brv.begin(); r != brv.end(); ++r) {
-      m_names.push_back(r->Qname());
-      fseq1_t *s;
-      
-      s = &m_seqs[n_seqs];
-      
-      s->seq   = strdup(r->Sequence().c_str());
-      s->qual  = strdup(r->Qualities().c_str());
+        s->seq  = strdup(r->Seq.c_str());
+        s->qual = strdup(r->Qual.c_str());
 
-      s->l_seq = r->Sequence().length();
-      size += m_seqs[n_seqs++].l_seq;
+        s->l_seq = r->Seq.length();
+        size += m_seqs[n_seqs++].l_seq;
+      }
+    }
+    void FermiAssembler::AddReads(const BamRecordVector &brv)
+    {
+      // alloc the memory
+      m_seqs =
+        (fseq1_t *)realloc(m_seqs, (n_seqs + brv.size()) * sizeof(fseq1_t));
+
+      uint64_t size = 0;
+      for (BamRecordVector::const_iterator r = brv.begin(); r != brv.end();
+           ++r) {
+        m_names.push_back(r->Qname());
+        fseq1_t *s;
+
+        s = &m_seqs[n_seqs];
+
+        s->seq  = strdup(r->Sequence().c_str());
+        s->qual = strdup(r->Qualities().c_str());
+
+        s->l_seq = r->Sequence().length();
+        size += m_seqs[n_seqs++].l_seq;
+      }
     }
-    
-  }
 
-  void FermiAssembler::ClearContigs() {
-    fml_utg_destroy(n_utg, m_utgs);  
-    m_utgs = 0;
-    n_utg = 0;
-  }
+    void FermiAssembler::ClearContigs()
+    {
+      fml_utg_destroy(n_utg, m_utgs);
+      m_utgs = 0;
+      n_utg  = 0;
+    }
 
-  void FermiAssembler::ClearReads() {  
-    if (!m_seqs)
-      return; //already cleared
-
-    for (size_t i = 0; i < n_seqs; ++i) {
-      fseq1_t * s = &m_seqs[i];
-      if (s->qual)
-       free(s->qual); 
-      s->qual = NULL;
-      if (s->seq)
-	free(s->seq);
-      s->seq = NULL;
+    void FermiAssembler::ClearReads()
+    {
+      if (!m_seqs) return; // already cleared
+
+      for (size_t i = 0; i < n_seqs; ++i) {
+        fseq1_t *s = &m_seqs[i];
+        if (s->qual) free(s->qual);
+        s->qual = NULL;
+        if (s->seq) free(s->seq);
+        s->seq = NULL;
+      }
+      free(m_seqs);
+      m_seqs = NULL;
     }
-    free(m_seqs);
-    m_seqs = NULL;
-      
-  }
 
-  void FermiAssembler::CorrectReads() {  
-    fml_correct(&opt, n_seqs, m_seqs);
-  }
+    void FermiAssembler::CorrectReads() { fml_correct(&opt, n_seqs, m_seqs); }
 
-  void FermiAssembler::CorrectAndFilterReads() {  
-    fml_fltuniq(&opt, n_seqs, m_seqs);
-  }
+    void FermiAssembler::CorrectAndFilterReads()
+    {
+      fml_fltuniq(&opt, n_seqs, m_seqs);
+    }
 
-  void FermiAssembler::PerformAssembly() {
-    m_utgs = fml_assemble(&opt, n_seqs, m_seqs, &n_utg); // assemble!
-  }
-  
-  std::vector<std::string> FermiAssembler::GetContigs() const {
-    std::vector<std::string> c;
-    for (size_t i = 0; i < n_utg; ++i)
-      c.push_back(std::string(m_utgs[i].seq));
-    return c;
+    void FermiAssembler::PerformAssembly()
+    {
+      m_utgs = fml_assemble(&opt, n_seqs, m_seqs, &n_utg); // assemble!
+    }
+
+    std::vector<std::string> FermiAssembler::GetContigs() const
+    {
+      std::vector<std::string> c;
+      for (size_t i = 0; i < n_utg; ++i)
+        c.push_back(std::string(m_utgs[i].seq));
+      return c;
+    }
+
+    /*void FermiAssembler::count() {
+
+      // initialize BFC options
+      uint64_t tot_len = 0;
+      for (int i = 0; i < n_seqs; ++i)
+        tot_len += m_seqs[i].l_seq; // compute total length
+      int l_pre = tot_len - 8 < 20? tot_len - 8 : 20;
+
+      //bfc_ch_t *ch = fml_count(n_seqs, m_seqs, opt.ec_k, 20, l_pre,
+      opt.n_threads);
+      //std::cerr << " ch->k " << ch->k << " ch->l_pre " << ch->l_pre <<
+      std::endl;
+
+      // directly from fml count
+      cnt_step_t cs;
+      cs.n_seqs = n_seqs, cs.seqs = m_seqs, cs.k = opt.ec_k, cs.q = 20;
+      cs.ch = bfc_ch_init(cs.k, l_pre);
+      }*/
+
+    UnalignedSequenceVector FermiAssembler::GetSequences() const
+    {
+      UnalignedSequenceVector r;
+      for (size_t i = 0; i < n_seqs; ++i) {
+        fseq1_t *s = &m_seqs[i];
+        UnalignedSequence read;
+        if (s->seq) read.Seq = (std::string(s->seq));
+        read.Name = m_names[i];
+        r.push_back(read);
+      }
+      return r;
+    }
   }
 
-  /*void FermiAssembler::count() {
-
-    // initialize BFC options
-    uint64_t tot_len = 0;
-    for (int i = 0; i < n_seqs; ++i) 
-      tot_len += m_seqs[i].l_seq; // compute total length
-    int l_pre = tot_len - 8 < 20? tot_len - 8 : 20;
-
-    //bfc_ch_t *ch = fml_count(n_seqs, m_seqs, opt.ec_k, 20, l_pre, opt.n_threads);
-    //std::cerr << " ch->k " << ch->k << " ch->l_pre " << ch->l_pre << std::endl;
-
-    // directly from fml count
-    cnt_step_t cs;
-    cs.n_seqs = n_seqs, cs.seqs = m_seqs, cs.k = opt.ec_k, cs.q = 20;
-    cs.ch = bfc_ch_init(cs.k, l_pre);
-    }*/
-
-  UnalignedSequenceVector FermiAssembler::GetSequences() const {
-    
-    UnalignedSequenceVector r;
-    for (size_t i = 0; i < n_seqs; ++i) {
-      fseq1_t * s = &m_seqs[i];
-      UnalignedSequence read;
-      if (s->seq)
-	read.Seq = (std::string(s->seq));
-      read.Name = m_names[i];
-      r.push_back(read);
+void
+SeqLib::FermiAssembler::WriteGFA(std::ostream &out)
+{
+  int i, j;
+  out << "H\tVN:Z:1.0" << std::endl;
+  for (i = 0; i < n_utg; ++i) {
+    const fml_utg_t *u = m_utgs + i;
+    out << "S\t" << i << "\t";
+    out << u->seq << "\tLN:i:" << u->len << "\tRC:i:" << u->nsr << "\tPD:Z:";
+    out << u->cov << std::endl;
+    for (j = 0; j < u->n_ovlp[0] + u->n_ovlp[1]; ++j) {
+      fml_ovlp_t *o = &u->ovlp[j];
+      if (i < o->id) {
+        out << "L\t" << i << "\t"
+            << "+-"[!o->from] << "\t" << o->id << "\t"
+            << "+-"[o->to] << "\t" << o->len << "M" << std::endl;
+      }
     }
-    return r;
   }
-  
 }


=====================================
src/GenomicRegion.cpp
=====================================
@@ -3,6 +3,9 @@
 #include <cassert>
 #include <stdexcept>
 #include <climits>
+#ifdef HAVE_C11
+#include <regex>
+#endif
 
 // 4 billion
 #define END_MAX 4000000000
@@ -61,9 +64,9 @@ int GenomicRegion::GetOverlap(const GenomicRegion& gr) const {
   }
 
   
-  std::string GenomicRegion::PointString() const {
+  std::string GenomicRegion::PointString(const BamHeader& h) const {
     std::stringstream out;
-    out << chrToString(chr) << ":" << SeqLib::AddCommas<int>(pos1) << "(" << strand << ")";
+    out << ChrName(h) << ":" << SeqLib::AddCommas<int>(pos1) << "(" << strand << ")";
     return out.str();
   }
 
@@ -112,11 +115,17 @@ bool GenomicRegion::operator>=(const GenomicRegion &b) const {
   return (*this > b || *this == b);
 }
 
-  std::string GenomicRegion::ToString() const {
+  /*  std::string GenomicRegion::ToString() const {
     return chrToString(chr) + ":" + SeqLib::AddCommas<int>(pos1) + "-" + AddCommas<int>(pos2) + "(" +
       strand + ")"; 
+      }*/
+
+  std::string GenomicRegion::ToString(const BamHeader& h) const {
+    return ChrName(h) + ":" + SeqLib::AddCommas<int>(pos1) + "-" + AddCommas<int>(pos2) + "(" +
+      strand + ")"; 
   }
 
+
 std::ostream& operator<<(std::ostream& out, const GenomicRegion& gr) {
   out << gr.chrToString(gr.chr) << ":" << SeqLib::AddCommas<int>(gr.pos1) << "-" << AddCommas<int>(gr.pos2) << "(" << 
     gr.strand << ")"; 
@@ -269,6 +278,16 @@ int32_t GenomicRegion::DistanceBetweenEnds(const GenomicRegion &gr) const {
       return;
     } else {
       chr = hdr.Name2ID(tchr); //bam_name2id(hdr.get(), tchr.c_str());
+
+      // if tchr was not found in sequence dictionary, it's possible that we're 
+      // specifying our contigs in b37 style whereas dict is hg** style;
+      // let's attempt to automatically convert [0-9XY]+ -> chr[0-9XY]+
+#ifdef HAVE_C11
+      static std::regex b37_regex("[0-9XY]+");
+      if(chr == -1 && std::regex_match(tchr, b37_regex)) {
+	 chr = hdr.Name2ID("chr" + tchr);
+      }
+#endif
     }
   }
 }


=====================================
src/Makefile.in
=====================================
@@ -1,7 +1,7 @@
-# Makefile.in generated by automake 1.15 from Makefile.am.
+# Makefile.in generated by automake 1.16.2 from Makefile.am.
 # @configure_input@
 
-# Copyright (C) 1994-2014 Free Software Foundation, Inc.
+# Copyright (C) 1994-2020 Free Software Foundation, Inc.
 
 # This Makefile.in is free software; the Free Software Foundation
 # gives unlimited permission to copy and/or distribute it,
@@ -86,7 +86,7 @@ POST_INSTALL = :
 NORMAL_UNINSTALL = :
 PRE_UNINSTALL = :
 POST_UNINSTALL = :
-subdir = src
+subdir = SeqLib/src
 ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
 am__aclocal_m4_deps = $(top_srcdir)/configure.ac
 am__configure_deps = $(am__aclocal_m4_deps) $(CONFIGURE_DEPENDENCIES) \
@@ -132,7 +132,22 @@ am__v_at_0 = @
 am__v_at_1 = 
 DEFAULT_INCLUDES = -I. at am__isrc@ -I$(top_builddir)
 depcomp = $(SHELL) $(top_srcdir)/depcomp
-am__depfiles_maybe = depfiles
+am__maybe_remake_depfiles = depfiles
+am__depfiles_remade = ./$(DEPDIR)/libseqlib_a-BFC.Po \
+	./$(DEPDIR)/libseqlib_a-BWAWrapper.Po \
+	./$(DEPDIR)/libseqlib_a-BamHeader.Po \
+	./$(DEPDIR)/libseqlib_a-BamReader.Po \
+	./$(DEPDIR)/libseqlib_a-BamRecord.Po \
+	./$(DEPDIR)/libseqlib_a-BamWriter.Po \
+	./$(DEPDIR)/libseqlib_a-FastqReader.Po \
+	./$(DEPDIR)/libseqlib_a-FermiAssembler.Po \
+	./$(DEPDIR)/libseqlib_a-GenomicRegion.Po \
+	./$(DEPDIR)/libseqlib_a-ReadFilter.Po \
+	./$(DEPDIR)/libseqlib_a-RefGenome.Po \
+	./$(DEPDIR)/libseqlib_a-SeqPlot.Po \
+	./$(DEPDIR)/libseqlib_a-jsoncpp.Po \
+	./$(DEPDIR)/libseqlib_a-ssw.Po \
+	./$(DEPDIR)/libseqlib_a-ssw_cpp.Po
 am__mv = mv -f
 AM_V_lt = $(am__v_lt_ at AM_V@)
 am__v_lt_ = $(am__v_lt_ at AM_DEFAULT_V@)
@@ -304,16 +319,16 @@ $(srcdir)/Makefile.in: @MAINTAINER_MODE_TRUE@ $(srcdir)/Makefile.am  $(am__confi
 	      exit 1;; \
 	  esac; \
 	done; \
-	echo ' cd $(top_srcdir) && $(AUTOMAKE) --foreign src/Makefile'; \
+	echo ' cd $(top_srcdir) && $(AUTOMAKE) --foreign SeqLib/src/Makefile'; \
 	$(am__cd) $(top_srcdir) && \
-	  $(AUTOMAKE) --foreign src/Makefile
+	  $(AUTOMAKE) --foreign SeqLib/src/Makefile
 Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status
 	@case '$?' in \
 	  *config.status*) \
 	    cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh;; \
 	  *) \
-	    echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe)'; \
-	    cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe);; \
+	    echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__maybe_remake_depfiles)'; \
+	    cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__maybe_remake_depfiles);; \
 	esac;
 
 $(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES)
@@ -339,21 +354,27 @@ mostlyclean-compile:
 distclean-compile:
 	-rm -f *.tab.c
 
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BFC.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BWAWrapper.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BamHeader.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BamReader.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BamRecord.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BamWriter.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-FastqReader.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-FermiAssembler.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-GenomicRegion.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-ReadFilter.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-RefGenome.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-SeqPlot.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-jsoncpp.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-ssw.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-ssw_cpp.Po at am__quote@
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BFC.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BWAWrapper.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BamHeader.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BamReader.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BamRecord.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-BamWriter.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-FastqReader.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-FermiAssembler.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-GenomicRegion.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-ReadFilter.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-RefGenome.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-SeqPlot.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-jsoncpp.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-ssw.Po at am__quote@ # am--include-marker
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/libseqlib_a-ssw_cpp.Po at am__quote@ # am--include-marker
+
+$(am__depfiles_remade):
+	@$(MKDIR_P) $(@D)
+	@echo '# dummy' >$@-t && $(am__mv) $@-t $@
+
+am--depfiles: $(am__depfiles_remade)
 
 .c.o:
 @am__fastdepCC_TRUE@	$(AM_V_CC)$(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $<
@@ -645,7 +666,10 @@ cscopelist-am: $(am__tagged_files)
 distclean-tags:
 	-rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags
 
-distdir: $(DISTFILES)
+distdir: $(BUILT_SOURCES)
+	$(MAKE) $(AM_MAKEFLAGS) distdir-am
+
+distdir-am: $(DISTFILES)
 	@srcdirstrip=`echo "$(srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \
 	topsrcdirstrip=`echo "$(top_srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \
 	list='$(DISTFILES)'; \
@@ -714,7 +738,21 @@ clean: clean-am
 clean-am: clean-generic clean-noinstLIBRARIES mostlyclean-am
 
 distclean: distclean-am
-	-rm -rf ./$(DEPDIR)
+		-rm -f ./$(DEPDIR)/libseqlib_a-BFC.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-BWAWrapper.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-BamHeader.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-BamReader.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-BamRecord.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-BamWriter.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-FastqReader.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-FermiAssembler.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-GenomicRegion.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-ReadFilter.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-RefGenome.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-SeqPlot.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-jsoncpp.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-ssw.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-ssw_cpp.Po
 	-rm -f Makefile
 distclean-am: clean-am distclean-compile distclean-generic \
 	distclean-tags
@@ -760,7 +798,21 @@ install-ps-am:
 installcheck-am:
 
 maintainer-clean: maintainer-clean-am
-	-rm -rf ./$(DEPDIR)
+		-rm -f ./$(DEPDIR)/libseqlib_a-BFC.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-BWAWrapper.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-BamHeader.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-BamReader.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-BamRecord.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-BamWriter.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-FastqReader.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-FermiAssembler.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-GenomicRegion.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-ReadFilter.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-RefGenome.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-SeqPlot.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-jsoncpp.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-ssw.Po
+	-rm -f ./$(DEPDIR)/libseqlib_a-ssw_cpp.Po
 	-rm -f Makefile
 maintainer-clean-am: distclean-am maintainer-clean-generic
 
@@ -780,18 +832,18 @@ uninstall-am:
 
 .MAKE: install-am install-strip
 
-.PHONY: CTAGS GTAGS TAGS all all-am check check-am clean clean-generic \
-	clean-noinstLIBRARIES cscopelist-am ctags ctags-am distclean \
-	distclean-compile distclean-generic distclean-tags distdir dvi \
-	dvi-am html html-am info info-am install install-am \
-	install-data install-data-am install-dvi install-dvi-am \
-	install-exec install-exec-am install-html install-html-am \
-	install-info install-info-am install-man install-pdf \
-	install-pdf-am install-ps install-ps-am install-strip \
-	installcheck installcheck-am installdirs maintainer-clean \
-	maintainer-clean-generic mostlyclean mostlyclean-compile \
-	mostlyclean-generic pdf pdf-am ps ps-am tags tags-am uninstall \
-	uninstall-am
+.PHONY: CTAGS GTAGS TAGS all all-am am--depfiles check check-am clean \
+	clean-generic clean-noinstLIBRARIES cscopelist-am ctags \
+	ctags-am distclean distclean-compile distclean-generic \
+	distclean-tags distdir dvi dvi-am html html-am info info-am \
+	install install-am install-data install-data-am install-dvi \
+	install-dvi-am install-exec install-exec-am install-html \
+	install-html-am install-info install-info-am install-man \
+	install-pdf install-pdf-am install-ps install-ps-am \
+	install-strip installcheck installcheck-am installdirs \
+	maintainer-clean maintainer-clean-generic mostlyclean \
+	mostlyclean-compile mostlyclean-generic pdf pdf-am ps ps-am \
+	tags tags-am uninstall uninstall-am
 
 .PRECIOUS: Makefile
 


=====================================
src/ReadFilter.cpp
=====================================
@@ -263,20 +263,27 @@ bool ReadFilter::isReadOverlappingRegion(const BamRecord &r) const {
       
     }
     
-    // check that there is at least one non-excluder region. 
-    // if not, give global includer
-    bool has_includer = false;
-    for (std::vector<ReadFilter>::const_iterator kk = m_regions.begin(); kk != m_regions.end(); ++kk) 
-      if (!kk->excluder)
-	has_includer = true;
-    if (!has_includer) {
-      ReadFilter mr;
-      mr.m_abstract_rules.push_back(rule_all);
-      mr.id = "WG_includer";
-      m_regions.push_back(mr);
-    }
-    
+    // make sure that there is at least one includer region
+    this->CheckHasIncluder();
+
   }
+
+    void ReadFilterCollection::CheckHasIncluder() {
+
+      // check that there is at least one non-excluder region. 
+      // if not, give global includer
+      bool has_includer = false;
+      for (std::vector<ReadFilter>::const_iterator kk = m_regions.begin(); kk != m_regions.end(); ++kk) 
+	if (!kk->excluder)
+	  has_includer = true;
+      if (!has_includer) {
+	ReadFilter mr;
+	mr.m_abstract_rules.push_back(rule_all);
+	mr.id = "WG_includer";
+	m_regions.push_back(mr);
+      }
+
+    }
   
   void ReadFilter::setRegions(const GRC& g) {
     m_grv = g;


=====================================
src/seqtools/Makefile.in
=====================================
@@ -1,7 +1,7 @@
-# Makefile.in generated by automake 1.15 from Makefile.am.
+# Makefile.in generated by automake 1.16.1 from Makefile.am.
 # @configure_input@
 
-# Copyright (C) 1994-2014 Free Software Foundation, Inc.
+# Copyright (C) 1994-2018 Free Software Foundation, Inc.
 
 # This Makefile.in is free software; the Free Software Foundation
 # gives unlimited permission to copy and/or distribute it,
@@ -119,7 +119,8 @@ am__v_at_0 = @
 am__v_at_1 = 
 DEFAULT_INCLUDES = -I. at am__isrc@ -I$(top_builddir)
 depcomp = $(SHELL) $(top_srcdir)/depcomp
-am__depfiles_maybe = depfiles
+am__maybe_remake_depfiles = depfiles
+am__depfiles_remade = ./$(DEPDIR)/seqtools-seqtools.Po
 am__mv = mv -f
 AM_V_lt = $(am__v_lt_ at AM_V@)
 am__v_lt_ = $(am__v_lt_ at AM_DEFAULT_V@)
@@ -284,8 +285,8 @@ Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status
 	  *config.status*) \
 	    cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh;; \
 	  *) \
-	    echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe)'; \
-	    cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe);; \
+	    echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__maybe_remake_depfiles)'; \
+	    cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__maybe_remake_depfiles);; \
 	esac;
 
 $(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES)
@@ -349,7 +350,13 @@ mostlyclean-compile:
 distclean-compile:
 	-rm -f *.tab.c
 
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/seqtools-seqtools.Po at am__quote@
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/seqtools-seqtools.Po at am__quote@ # am--include-marker
+
+$(am__depfiles_remade):
+	@$(MKDIR_P) $(@D)
+	@echo '# dummy' >$@-t && $(am__mv) $@-t $@
+
+am--depfiles: $(am__depfiles_remade)
 
 .cpp.o:
 @am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXXCOMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $<
@@ -431,7 +438,10 @@ cscopelist-am: $(am__tagged_files)
 distclean-tags:
 	-rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags
 
-distdir: $(DISTFILES)
+distdir: $(BUILT_SOURCES)
+	$(MAKE) $(AM_MAKEFLAGS) distdir-am
+
+distdir-am: $(DISTFILES)
 	@srcdirstrip=`echo "$(srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \
 	topsrcdirstrip=`echo "$(top_srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \
 	list='$(DISTFILES)'; \
@@ -503,7 +513,7 @@ clean: clean-am
 clean-am: clean-binPROGRAMS clean-generic mostlyclean-am
 
 distclean: distclean-am
-	-rm -rf ./$(DEPDIR)
+		-rm -f ./$(DEPDIR)/seqtools-seqtools.Po
 	-rm -f Makefile
 distclean-am: clean-am distclean-compile distclean-generic \
 	distclean-tags
@@ -549,7 +559,7 @@ install-ps-am:
 installcheck-am:
 
 maintainer-clean: maintainer-clean-am
-	-rm -rf ./$(DEPDIR)
+		-rm -f ./$(DEPDIR)/seqtools-seqtools.Po
 	-rm -f Makefile
 maintainer-clean-am: distclean-am maintainer-clean-generic
 
@@ -569,7 +579,7 @@ uninstall-am: uninstall-binPROGRAMS
 
 .MAKE: install-am install-strip
 
-.PHONY: CTAGS GTAGS TAGS all all-am check check-am clean \
+.PHONY: CTAGS GTAGS TAGS all all-am am--depfiles check check-am clean \
 	clean-binPROGRAMS clean-generic cscopelist-am ctags ctags-am \
 	distclean distclean-compile distclean-generic distclean-tags \
 	distdir dvi dvi-am html html-am info info-am install \



View it on GitLab: https://salsa.debian.org/med-team/libseqlib/-/commit/438a32ffba00bda61a04d0ae06dad7fa101c0f37

-- 
View it on GitLab: https://salsa.debian.org/med-team/libseqlib/-/commit/438a32ffba00bda61a04d0ae06dad7fa101c0f37
You're receiving this email because of your account on salsa.debian.org.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20201217/2e20ffcc/attachment-0001.html>


More information about the debian-med-commit mailing list