Merge branch 'master' of github.com:walaj/SeqLib into jh_nov2018

Author: julianhess

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

@@ -11,10 +11,11 @@ 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
@@ -30,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; }

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; }
 
@@ -636,9 +639,9 @@ class BamRecord {
     uint32_t* c = bam_get_cigar(b);
     int32_t p = 0;
     for (size_t i = 0; i < b->core.n_cigar; ++i) {
-      if ( (bam_cigar_opchr(c[i]) == 'S') || (bam_cigar_opchr(c[i]) == 'H'))
+      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;
@@ -826,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/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/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/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/UnalignedSequence.h

@@ -16,6 +16,7 @@ extern "C" {
 
 #include <cstring>
 #include <vector>
+#include <iostream>
 
 namespace SeqLib {
 
@@ -30,27 +31,39 @@ 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 
      * @param s Sequence, stored as ACTG or N characters
      * @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