@@ -10,7 +10,7 @@ Description: GSNAP and GMAP are a pair of tools to align short-read

     to work with GMAP and GSNAP from within R. In addition, it provides 

     methods to tally alignment results on a per-nucleotide basis using 

     the bam_tally tool.

-Version: 1.3.9

+Version: 1.3.10

 Depends: R (>= 2.15.0), methods, GenomicRanges

 Imports: IRanges, Rsamtools (>= 1.7.4), rtracklayer (>= 1.17.15), GenomicRanges,

          GenomicFeatures, Biostrings, VariantAnnotation (>= 1.7.34), tools,

@@ -8,6 +8,7 @@

 setClass("BamTallyParam",

          representation(genome = "GmapGenome",

                         which = "GenomicRanges",

+                        desired_read_group = "characterORNULL",

                         minimum_mapq = "integer",

                         concordant_only = "logical",

                         unique_only = "logical",

@@ -32,13 +33,15 @@ normArgWhich <- function(x, genome) {

 }

 BamTallyParam <- function(genome, which = GRanges(),

-                          minimum_mapq = 0L,

+                          desired_read_group = NULL, minimum_mapq = 0L,

                           concordant_only = FALSE, unique_only = FALSE,

                           primary_only = FALSE, ignore_duplicates = FALSE,

                           min_depth = 0L, variant_strand = 0L,

                           ignore_query_Ns = FALSE,

                           indels = FALSE)

 {

+  if (!is.null(desired_read_group) && !isSingleString(desired_read_group))

+    stop("'desired_read_group' must be NULL or a single, non-NA string")

   if (!isSingleNumber(minimum_mapq) || minimum_mapq < 0)

     stop("minimum_mapq must be a single, non-negative, non-NA number")

   if (!isTRUEorFALSE(concordant_only))

@@ -120,11 +120,11 @@ normalizeIndelAlleles <- function(x, genome) {

   is.indel <- nchar(ref(x)) == 0L | nchar(alt(x)) == 0L

   if (any(is.indel)) {

     indels <- x[is.indel]

-    indels <- flank(indels, 1)

-    anchor <- getSeq(genome, indels)

-    ref(indels) <- paste0(anchor, ref(indels))

-    alt(indels) <- paste0(anchor, alt(indels))

-    x[is.indel] <- indels

+    flanks <- flank(indels, 1)

+    anchor <- getSeq(genome, flanks)

+    ref(x)[is.indel] <- paste0(anchor, ref(indels))

+    alt(x)[is.indel] <- paste0(anchor, alt(indels))

+    ranges(x)[is.indel] <- resize(ranges(indels), nchar(ref(x)[is.indel]), "end")

   }

   x

 }

@@ -136,7 +136,7 @@ normalizeIndelAlleles <- function(x, genome) {

 normArgSingleInteger <- function(x) {

   name <- deparse(substitute(x))

   x <- as.integer(x)

-  if (!IRanges:::isSingleInteger(x))

+  if (!isSingleInteger(x))

     stop("'", name, "' should be a single, non-NA integer")

   x

 }

@@ -149,7 +149,8 @@ normArgTRUEorFALSE <- function(x) {

 .bam_tally_C <- function(bamreader, genome_dir = NULL, db = NULL,

                          which = NULL, read_pos_breaks = NULL,

-                         high_base_quality = 0L, alloclength = 200000L,

+                         high_base_quality = 0L, desired_read_group = NULL,

+                         alloclength = 200000L,

                          minimum_mapq = 0L, good_unique_mapq = 35L,

                          maximum_nhits = 1000000L,

                          concordant_only = FALSE, unique_only = FALSE,

@@ -164,10 +165,12 @@ normArgTRUEorFALSE <- function(x) {

   if (length(which) > 0L) {

     which <- list(as.character(seqnames(which)), start(which), end(which))

   } else which <- NULL

-  if (!is.null(genome_dir) && !IRanges:::isSingleString(genome_dir))

+  if (!is.null(genome_dir) && !isSingleString(genome_dir))

     stop("'genome_dir' must be NULL or a single, non-NA string")

-  if (!is.null(db) && !IRanges:::isSingleString(db))

+  if (!is.null(db) && !isSingleString(db))

     stop("'db' must be NULL or a single, non-NA string")

+  if (!is.null(desired_read_group) && !isSingleString(desired_read_group))

+    stop("'desired_read_group' must be NULL or a single, non-NA string")

   if (!is.null(read_pos_breaks)) {

     read_pos_breaks <- as.integer(read_pos_breaks)

     if (any(is.na(read_pos_breaks)))

@@ -178,6 +181,7 @@ normArgTRUEorFALSE <- function(x) {

       stop("'read_pos_breaks' must be sorted")

   }

   .Call(R_Bamtally_iit, bamreader@.extptr, genome_dir, db, which,

+        desired_read_group,

         normArgSingleInteger(alloclength),

         normArgSingleInteger(minimum_mapq),

         normArgSingleInteger(good_unique_mapq),

@@ -14,14 +14,16 @@

 }

 \usage{

-  BamTallyParam(genome, which = RangesList(), read_pos_breaks = NULL,

+  BamTallyParam(genome, which = RangesList(),

+                desired_read_group = NULL,

+                read_pos_breaks = NULL,

                 high_base_quality = 0L,

                 minimum_mapq = 0L,

                 concordant_only = FALSE, unique_only = FALSE,

                 primary_only = FALSE, ignore_duplicates = FALSE,

                 min_depth = 0L, variant_strand = 0L,

                 ignore_query_Ns = FALSE,

-                indels = FALSE)

+                indels = TRUE)

 }

 \arguments{

   \item{genome}{A \code{GmapGenome} object, or something coercible to one.}

@@ -29,6 +31,8 @@

     one that limits the tally to that range or set of ranges. By

     default, the entire genome is processed.

   }

+  \item{desired_read_group}{The name of the read group to which to limit

+    the tallying; if not NULL, must be a single, non-NA string.}

   \item{read_pos_breaks}{The breaks, like those passed to \code{\link{cut}}

     for aggregating the per-read position counts. If \code{NULL}, no per-cycle

     counts are returned.}

@@ -54,7 +58,12 @@

     is a good way to save resources.}

   \item{ignore_query_Ns}{Whether to ignore the N base pairs when

     counting. Can save a lot of resources when processing low quality data.}

-  \item{indels}{Whether to return indel counts; not supported yet.}

+  \item{indels}{Whether to return indel counts. The \code{ref} and

+    \code{alt} columns in the returned \code{VRanges} conform to VCF

+    conventions; i.e., the first base upstream is included. The range

+    always spans the sequence in \code{ref}; so e.g. a deletion extends

+    one nt upstream of the actual deleted sequence.

+  }

 }

 \seealso{

@@ -49,7 +49,7 @@ gstruct/Makefile: gstruct/configure ${R_SRC_DIR}/samtools

 	            --prefix=${PREFIX} --includedir=${GSTRUCT_INCLUDE_DIR} \

 	            --libdir=${PREFIX}/${LIBnn} \

 	            --with-samtools=${R_SRC_DIR}/samtools \

-		    --disable-binaries --disable-maintainer-mode #\

+		    --disable-binaries #--disable-maintainer-mode #\

 ## does not appear to be a true dependency yet

 ##	            --with-gmap=${PREFIX}/bin

@@ -7,7 +7,7 @@

 static const R_CallMethodDef callMethods[] = {

   /* bamtally.c */

-  CALLMETHOD_DEF(R_Bamtally_iit, 18),

+  CALLMETHOD_DEF(R_Bamtally_iit, 19),

   CALLMETHOD_DEF(R_tally_iit_parse, 4),

   /* bamreader.c */

@@ -13,7 +13,7 @@

 SEXP

 R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,

-                SEXP which_R,

+                SEXP which_R, SEXP desired_read_group_R,

                 SEXP alloclength_R,

                 SEXP minimum_mapq_R, SEXP good_unique_mapq_R,

                 SEXP maximum_nhits_R,

@@ -30,6 +30,9 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,

     genome_dir_R == R_NilValue ? NULL : CHAR(asChar(genome_dir_R));

   const char *db = CHAR(asChar(db_R));

   int alloclength = asInteger(alloclength_R);

+  const char *desired_read_group =

+    desired_read_group_R == R_NilValue ? NULL :

+    CHAR(asChar(desired_read_group_R));

   int minimum_mapq = asInteger(minimum_mapq_R);

   int good_unique_mapq = asInteger(good_unique_mapq_R);

   int maximum_nhits = asInteger(maximum_nhits_R);

@@ -57,9 +60,12 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,

     end = asInteger(VECTOR_ELT(which_R, 2));

   }

-  IIT_T tally_iit = Bamtally_iit(bamreader, (char *)chr, start, end,

-                                 genome, chromosome_iit,

-                                 alloclength, minimum_mapq, good_unique_mapq,

+  IIT_T tally_iit = Bamtally_iit(bamreader, (char *)chr, 

+                                 /* TODO: bam_lacks_chr */ NULL,

+                                 start, end,

+                                 genome, chromosome_iit, alloclength,

+                                 (char *)desired_read_group,

+                                 minimum_mapq, good_unique_mapq,

                                  maximum_nhits, need_concordant_p,

                                  need_unique_p, need_primary_p,

                                  ignore_duplicates_p,

@@ -10,7 +10,7 @@ R_Bamread_new (SEXP bamfile_R);

 SEXP

 R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,

-                SEXP which_R,

+                SEXP which_R, SEXP desired_read_group_R,

                 SEXP alloclength_R,

                 SEXP minimum_mapq_R, SEXP good_unique_mapq_R,

                 SEXP maximum_nhits_R,

@@ -51,13 +51,13 @@ build_triplet = @build@

 host_triplet = @host@

 target_triplet = @target@

 subdir = .

-DIST_COMMON = INSTALL NEWS README AUTHORS ChangeLog \

+DIST_COMMON = INSTALL NEWS README AUTHORS ChangeLog COPYING \

 	$(srcdir)/Makefile.in $(srcdir)/Makefile.am \

 	$(top_srcdir)/configure $(am__configure_deps) \

-	$(top_srcdir)/lib/gstruct.pc.in COPYING config/compile \

-	config/config.guess config/config.sub config/depcomp \

-	config/install-sh config/missing config/ltmain.sh \

-	$(top_srcdir)/config/compile $(top_srcdir)/config/config.guess \

+	$(top_srcdir)/lib/gstruct.pc.in config/compile \

+	config/config.guess config/config.sub config/install-sh \

+	config/missing config/ltmain.sh $(top_srcdir)/config/compile \

+	$(top_srcdir)/config/config.guess \

 	$(top_srcdir)/config/config.sub \

 	$(top_srcdir)/config/install-sh $(top_srcdir)/config/ltmain.sh \

 	$(top_srcdir)/config/missing

@@ -231,7 +231,6 @@ LIBTOOL = @LIBTOOL@

 LIPO = @LIPO@

 LN_S = @LN_S@

 LTLIBOBJS = @LTLIBOBJS@

-MAINT = @MAINT@

 MAKEINFO = @MAKEINFO@

 MANIFEST_TOOL = @MANIFEST_TOOL@

 MKDIR_P = @MKDIR_P@

@@ -329,7 +328,7 @@ all: all-recursive

 .SUFFIXES:

 am--refresh: Makefile

 	@:

-$(srcdir)/Makefile.in: @MAINTAINER_MODE_TRUE@ $(srcdir)/Makefile.am  $(am__configure_deps)

+$(srcdir)/Makefile.in:  $(srcdir)/Makefile.am  $(am__configure_deps)

 	@for dep in $?; do \

 	  case '$(am__configure_deps)' in \

 	    *$$dep*) \

@@ -356,9 +355,9 @@ Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status

 $(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES)

 	$(SHELL) ./config.status --recheck

-$(top_srcdir)/configure: @MAINTAINER_MODE_TRUE@ $(am__configure_deps)

+$(top_srcdir)/configure:  $(am__configure_deps)

 	$(am__cd) $(srcdir) && $(AUTOCONF)

-$(ACLOCAL_M4): @MAINTAINER_MODE_TRUE@ $(am__aclocal_m4_deps)

+$(ACLOCAL_M4):  $(am__aclocal_m4_deps)

 	$(am__cd) $(srcdir) && $(ACLOCAL) $(ACLOCAL_AMFLAGS)

 $(am__aclocal_m4_deps):

 lib/gstruct-${LIBGSTRUCT_API_VERSION}.pc: $(top_builddir)/config.status $(top_srcdir)/lib/gstruct.pc.in

@@ -1,22 +1,2 @@

-

-Usage:

-

-bam_singlemap -d <genome> <bam> | samtools view -Sb - > <new bam>

-

-bam_tally -d <genome> <bam>

-bam_gstruct -d <genome> <bam> | iit_store -o <pred>

-bam_print -d <genome> -m <genes>,<pred> <bam> <coords>

-

-

-genecompare -d <genome> -m <genes> <pred.iit> > x, or

-genecompare <goldstd.iit> <pred.iit> > x

-

-

-Depends on samtools and gmap packages.

-

-gfusion depends on gmap and qsub

-

-

-

@@ -1 +1 @@

-2013-02-14

\ No newline at end of file

+2013-09-12

\ No newline at end of file

@@ -596,42 +596,6 @@ fi

 rmdir .tst 2>/dev/null

 AC_SUBST([am__leading_dot])])

-# Add --enable-maintainer-mode option to configure.         -*- Autoconf -*-

-# From Jim Meyering

-

-# Copyright (C) 1996-2013 Free Software Foundation, Inc.

-#

-# This file is free software; the Free Software Foundation

-# gives unlimited permission to copy and/or distribute it,

-# with or without modifications, as long as this notice is preserved.

-

-# AM_MAINTAINER_MODE([DEFAULT-MODE])

-# ----------------------------------

-# Control maintainer-specific portions of Makefiles.

-# Default is to disable them, unless 'enable' is passed literally.

-# For symmetry, 'disable' may be passed as well.  Anyway, the user

-# can override the default with the --enable/--disable switch.

-AC_DEFUN([AM_MAINTAINER_MODE],

-[m4_case(m4_default([$1], [disable]),

-       [enable], [m4_define([am_maintainer_other], [disable])],

-       [disable], [m4_define([am_maintainer_other], [enable])],

-       [m4_define([am_maintainer_other], [enable])

-        m4_warn([syntax], [unexpected argument to AM@&t@_MAINTAINER_MODE: $1])])

-AC_MSG_CHECKING([whether to enable maintainer-specific portions of Makefiles])

-  dnl maintainer-mode's default is 'disable' unless 'enable' is passed

-  AC_ARG_ENABLE([maintainer-mode],

-    [AS_HELP_STRING([--]am_maintainer_other[-maintainer-mode],

-      am_maintainer_other[ make rules and dependencies not useful

-      (and sometimes confusing) to the casual installer])],

-    [USE_MAINTAINER_MODE=$enableval],

-    [USE_MAINTAINER_MODE=]m4_if(am_maintainer_other, [enable], [no], [yes]))

-  AC_MSG_RESULT([$USE_MAINTAINER_MODE])

-  AM_CONDITIONAL([MAINTAINER_MODE], [test $USE_MAINTAINER_MODE = yes])

-  MAINT=$MAINTAINER_MODE_TRUE

-  AC_SUBST([MAINT])dnl

-]

-)

-

 # Check to see how 'make' treats includes.	            -*- Autoconf -*-

 # Copyright (C) 2001-2013 Free Software Foundation, Inc.

@@ -1,6 +1,6 @@

 #! /bin/sh

 # Guess values for system-dependent variables and create Makefiles.

-# Generated by GNU Autoconf 2.69 for gstruct 2013-02-14.

+# Generated by GNU Autoconf 2.69 for gstruct 2013-09-12.

 #

 # Report bugs to <Thomas Wu <twu@gene.com>>.

 #

@@ -590,8 +590,8 @@ MAKEFLAGS=

 # Identity of this package.

 PACKAGE_NAME='gstruct'

 PACKAGE_TARNAME='gstruct'

-PACKAGE_VERSION='2013-02-14'

-PACKAGE_STRING='gstruct 2013-02-14'

+PACKAGE_VERSION='2013-09-12'

+PACKAGE_STRING='gstruct 2013-09-12'

 PACKAGE_BUGREPORT='Thomas Wu <twu@gene.com>'

 PACKAGE_URL=''

@@ -684,9 +684,6 @@ MAINTAINER_FALSE

 MAINTAINER_TRUE

 FULLDIST_FALSE

 FULLDIST_TRUE

-MAINT

-MAINTAINER_MODE_FALSE

-MAINTAINER_MODE_TRUE

 AM_BACKSLASH

 AM_DEFAULT_VERBOSITY

 AM_DEFAULT_V

@@ -789,7 +786,6 @@ enable_option_checking

 enable_largefile

 enable_dependency_tracking

 enable_silent_rules

-enable_maintainer_mode

 enable_fulldist

 enable_maintainer

 enable_static_linking

@@ -1355,7 +1351,7 @@ if test "$ac_init_help" = "long"; then

   # Omit some internal or obsolete options to make the list less imposing.

   # This message is too long to be a string in the A/UX 3.1 sh.

   cat <<_ACEOF

-\`configure' configures gstruct 2013-02-14 to adapt to many kinds of systems.

+\`configure' configures gstruct 2013-09-12 to adapt to many kinds of systems.

 Usage: $0 [OPTION]... [VAR=VALUE]...

@@ -1426,7 +1422,7 @@ fi

 if test -n "$ac_init_help"; then

   case $ac_init_help in

-     short | recursive ) echo "Configuration of gstruct 2013-02-14:";;

+     short | recursive ) echo "Configuration of gstruct 2013-09-12:";;

    esac

   cat <<\_ACEOF

@@ -1441,9 +1437,6 @@ Optional Features:

                           speeds up one-time build

   --enable-silent-rules   less verbose build output (undo: "make V=1")

   --disable-silent-rules  verbose build output (undo: "make V=0")

-  --disable-maintainer-mode

-                          disable make rules and dependencies not useful (and

-                          sometimes confusing) to the casual installer

   --enable-fulldist       For use by program maintainer

   --enable-maintainer     For use by program maintainer

   --enable-static-linking Link binaries statically (default=no)

@@ -1545,7 +1538,7 @@ fi

 test -n "$ac_init_help" && exit $ac_status

 if $ac_init_version; then

   cat <<\_ACEOF

-gstruct configure 2013-02-14

+gstruct configure 2013-09-12

 generated by GNU Autoconf 2.69

 Copyright (C) 2012 Free Software Foundation, Inc.

@@ -2151,7 +2144,7 @@ cat >config.log <<_ACEOF

 This file contains any messages produced by compilers while

 running configure, to aid debugging if configure makes a mistake.

-It was created by gstruct $as_me 2013-02-14, which was

+It was created by gstruct $as_me 2013-09-12, which was

 generated by GNU Autoconf 2.69.  Invocation command line was

   $ $0 $@

@@ -2501,8 +2494,8 @@ ac_compiler_gnu=$ac_cv_c_compiler_gnu

 { $as_echo "$as_me:${as_lineno-$LINENO}: checking package version" >&5

 $as_echo_n "checking package version... " >&6; }

-{ $as_echo "$as_me:${as_lineno-$LINENO}: result: 2013-02-14" >&5

-$as_echo "2013-02-14" >&6; }

+{ $as_echo "$as_me:${as_lineno-$LINENO}: result: 2013-09-12" >&5

+$as_echo "2013-09-12" >&6; }

 ### Read defaults

 { $as_echo "$as_me:${as_lineno-$LINENO}: checking CONFIG_SITE" >&5

@@ -4256,8 +4249,9 @@ fi

 # Define the identity of the package.

- PACKAGE='gstruct'

- VERSION='2013-02-14'

+

+ PACKAGE=gstruct

+ VERSION=2013-09-12

 cat >>confdefs.h <<_ACEOF

@@ -4432,29 +4426,6 @@ fi

-{ $as_echo "$as_me:${as_lineno-$LINENO}: checking whether to enable maintainer-specific portions of Makefiles" >&5

-$as_echo_n "checking whether to enable maintainer-specific portions of Makefiles... " >&6; }

-    # Check whether --enable-maintainer-mode was given.

-if test "${enable_maintainer_mode+set}" = set; then :

-  enableval=$enable_maintainer_mode; USE_MAINTAINER_MODE=$enableval

-else

-  USE_MAINTAINER_MODE=yes

-fi

-

-  { $as_echo "$as_me:${as_lineno-$LINENO}: result: $USE_MAINTAINER_MODE" >&5

-$as_echo "$USE_MAINTAINER_MODE" >&6; }

-   if test $USE_MAINTAINER_MODE = yes; then

-  MAINTAINER_MODE_TRUE=

-  MAINTAINER_MODE_FALSE='#'

-else

-  MAINTAINER_MODE_TRUE='#'

-  MAINTAINER_MODE_FALSE=

-fi

-

-  MAINT=$MAINTAINER_MODE_TRUE

-

-

-

  if test "x$enable_fulldist" = xyes; then

   FULLDIST_TRUE=

   FULLDIST_FALSE='#'

@@ -15457,10 +15428,6 @@ else

   am__EXEEXT_FALSE=

 fi

-if test -z "${MAINTAINER_MODE_TRUE}" && test -z "${MAINTAINER_MODE_FALSE}"; then

-  as_fn_error $? "conditional \"MAINTAINER_MODE\" was never defined.

-Usually this means the macro was only invoked conditionally." "$LINENO" 5

-fi

 if test -z "${FULLDIST_TRUE}" && test -z "${FULLDIST_FALSE}"; then

   as_fn_error $? "conditional \"FULLDIST\" was never defined.

 Usually this means the macro was only invoked conditionally." "$LINENO" 5

@@ -15875,7 +15842,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1

 # report actual input values of CONFIG_FILES etc. instead of their

 # values after options handling.

 ac_log="

-This file was extended by gstruct $as_me 2013-02-14, which was

+This file was extended by gstruct $as_me 2013-09-12, which was

 generated by GNU Autoconf 2.69.  Invocation command line was

   CONFIG_FILES    = $CONFIG_FILES

@@ -15941,7 +15908,7 @@ _ACEOF

 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1

 ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"

 ac_cs_version="\\

-gstruct config.status 2013-02-14

+gstruct config.status 2013-09-12

 configured by $0, generated by GNU Autoconf 2.69,

   with options \\"\$ac_cs_config\\"

@@ -63,8 +63,7 @@ AC_SYS_LARGEFILE

 AC_ARG_PROGRAM

 #AM_INIT_AUTOMAKE([no-dependencies])

-AM_INIT_AUTOMAKE([1.10.0])

-AM_MAINTAINER_MODE([enable])

+AM_INIT_AUTOMAKE(AC_PACKAGE_NAME, AC_PACKAGE_VERSION)

 AM_CONDITIONAL(FULLDIST,test "x$enable_fulldist" = xyes)

 AC_ARG_ENABLE([fulldist],

@@ -360,7 +360,6 @@ LIBTOOL = @LIBTOOL@

 LIPO = @LIPO@

 LN_S = @LN_S@

 LTLIBOBJS = @LTLIBOBJS@

-MAINT = @MAINT@

 MAKEINFO = @MAKEINFO@

 MANIFEST_TOOL = @MANIFEST_TOOL@

 MKDIR_P = @MKDIR_P@

@@ -705,7 +704,7 @@ all: config.h

 .SUFFIXES:

 .SUFFIXES: .c .lo .o .obj

-$(srcdir)/Makefile.in: @MAINTAINER_MODE_TRUE@ $(srcdir)/Makefile.am  $(am__configure_deps)

+$(srcdir)/Makefile.in:  $(srcdir)/Makefile.am  $(am__configure_deps)

 	@for dep in $?; do \

 	  case '$(am__configure_deps)' in \

 	    *$$dep*) \

@@ -730,9 +729,9 @@ Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status

 $(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES)

 	cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh

-$(top_srcdir)/configure: @MAINTAINER_MODE_TRUE@ $(am__configure_deps)

+$(top_srcdir)/configure:  $(am__configure_deps)

 	cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh

-$(ACLOCAL_M4): @MAINTAINER_MODE_TRUE@ $(am__aclocal_m4_deps)

+$(ACLOCAL_M4):  $(am__aclocal_m4_deps)

 	cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh

 $(am__aclocal_m4_deps):

@@ -743,7 +742,7 @@ config.h: stamp-h1

 stamp-h1: $(srcdir)/config.h.in $(top_builddir)/config.status

 	@rm -f stamp-h1

 	cd $(top_builddir) && $(SHELL) ./config.status src/config.h

-$(srcdir)/config.h.in: @MAINTAINER_MODE_TRUE@ $(am__configure_deps) 

+$(srcdir)/config.h.in:  $(am__configure_deps) 

 	($(am__cd) $(top_srcdir) && $(AUTOHEADER))

 	rm -f stamp-h1

 	touch $@

@@ -1,4 +1,4 @@

-static char rcsid[] = "$Id: bam_fastq.c 87713 2013-03-01 18:32:34Z twu $";

+static char rcsid[] = "$Id: bam_fastq.c 108654 2013-09-19 23:11:00Z twu $";

 #ifdef HAVE_CONFIG_H

 #include <config.h>

 #endif

@@ -57,6 +57,7 @@ static char rcsid[] = "$Id: bam_fastq.c 87713 2013-03-01 18:32:34Z twu $";

 static char *output_root = "out";

 static char *error_root = NULL;

 static bool fastap = false;

+static bool add_bounds_p = false;

 static bool allp = true;

 static bool single_end_output_p = false; /* print paired-end BAM as single-end file */

 static int terminal_minlength = -1;

@@ -75,6 +76,7 @@ static struct option long_options[] = {

   {"output", required_argument, 0, 'o'},	/* output_root */

   {"error", required_argument, 0, 'e'},	/* error_root */

   {"fasta", no_argument, 0, 0},			/* fastap */

+  {"add-bounds", no_argument, 0, 0},		/* add_bounds_p */

   {"all", no_argument, 0, 'A'}, /* allp */

   {"single-end-output", no_argument, 0, '1'}, /* single_end_output_p */

   {"spliced", no_argument, 0, 0}, /* splicedp */

@@ -109,6 +111,7 @@ Usage: bam_fastq [options] [-o <output>] <bam file...>\n\

 \n\

 Options:\n\

   --fasta                         Writes to stdout in FASTA format\n\

+  --add-lh                        Add chromosomal low/high designation for concordant paired-end alignments\n\

   -o, --output=STRING             Writes FASTQ to <output>_1.fq and <output>_2.fq (default=out)\n\

   -1, --single-end-output         Prints reads as single-end (even if BAM files have paired-end reads)\n\

   -S, --terminal-minlength=INT    Print only reads with terminal ends (cigar S length) equal to\n\

@@ -197,11 +200,35 @@ struct T {

   char *queryseq;

   char *quality_string;

   bool splicedp;

+

+  char querystart_strand;

+  char *querystart_chr;

+  Genomicpos_T querystart_chrpos;

+

+  char querybreak1_strand;

+  char *querybreak1_chr;

+  Genomicpos_T querybreak1_chrpos;

+

+  char querybreak2_strand;

+  char *querybreak2_chr;

+  Genomicpos_T querybreak2_chrpos;

+

+  char queryend_strand;

+  char *queryend_chr;

+  Genomicpos_T queryend_chrpos;

 };

 static void

 Read_free (T *old) {

+  if ((*old)->querybreak1_chrpos != 0) {

+    FREE((*old)->querybreak1_chr);

+  }

+  if ((*old)->querybreak2_chrpos != 0) {

+    FREE((*old)->querybreak2_chr);

+  }

+  FREE((*old)->queryend_chr);

+  FREE((*old)->querystart_chr);

   FREE((*old)->queryseq);

   if ((*old)->quality_string != NULL) {

     FREE((*old)->quality_string);

@@ -211,6 +238,102 @@ Read_free (T *old) {

 }

+static void

+Read_add_chrinfo_querystart (Read_T read, int flag, char *chr, Genomicpos_T chrpos_low, Genomicpos_T chrpos_high) {

+

+  if (chr == NULL) {

+    read->querystart_strand = '*';

+    read->querystart_chr = (char *) CALLOC(1,sizeof(char));

+    read->querystart_chr[0] = '\0';

+    read->querystart_chrpos = 0;

+  } else {

+    read->querystart_chr = (char *) CALLOC(strlen(chr)+1,sizeof(char));

+    strcpy(read->querystart_chr,chr);

+    if (flag & QUERY_MINUSP) {

+      read->querystart_strand = '-';

+      read->querystart_chrpos = chrpos_high;

+    } else {

+      read->querystart_strand = '+';

+      read->querystart_chrpos = chrpos_low;

+    }

+  }

+

+  return;

+}

+

+static void

+Read_add_chrinfo_querybreak2 (Read_T read, int flag, char *chr, Genomicpos_T chrpos_low, Genomicpos_T chrpos_high) {

+

+  if (chr == NULL) {

+    abort();

+    read->querybreak2_strand = '*';

+    read->querybreak2_chr = (char *) CALLOC(1,sizeof(char));

+    read->querybreak2_chr[0] = '\0';

+    read->querybreak2_chrpos = 0;

+  } else {

+    read->querybreak2_chr = (char *) CALLOC(strlen(chr)+1,sizeof(char));

+    strcpy(read->querybreak2_chr,chr);

+    if (flag & QUERY_MINUSP) {

+      read->querybreak2_strand = '-';

+      read->querybreak2_chrpos = chrpos_high;

+    } else {

+      read->querybreak2_strand = '+';

+      read->querybreak2_chrpos = chrpos_low;

+    }

+  }

+

+  return;

+}

+

+

+static void

+Read_add_chrinfo_queryend (Read_T read, int flag, char *chr, Genomicpos_T chrpos_low, Genomicpos_T chrpos_high) {

+

+  if (chr == NULL) {

+    read->queryend_strand = '*';

+    read->queryend_chr = (char *) CALLOC(1,sizeof(char));

+    read->queryend_chr[0] = '\0';

+    read->queryend_chrpos = 0;

+  } else {

+    read->queryend_chr = (char *) CALLOC(strlen(chr)+1,sizeof(char));

+    strcpy(read->queryend_chr,chr);

+    if (flag & QUERY_MINUSP) {

+      read->queryend_strand = '-';

+      read->queryend_chrpos = chrpos_low;

+    } else {

+      read->queryend_strand = '+';

+      read->queryend_chrpos = chrpos_high;

+    }

+  }

+

+  return;

+}

+

+static void

+Read_add_chrinfo_querybreak1 (Read_T read, int flag, char *chr, Genomicpos_T chrpos_low, Genomicpos_T chrpos_high) {

+

+  if (chr == NULL) {

+    abort();

+    read->querybreak1_strand = '*';

+    read->querybreak1_chr = (char *) CALLOC(1,sizeof(char));

+    read->querybreak1_chr[0] = '\0';

+    read->querybreak1_chrpos = 0;

+  } else {

+    read->querybreak1_chr = (char *) CALLOC(strlen(chr)+1,sizeof(char));

+    strcpy(read->querybreak1_chr,chr);

+    if (flag & QUERY_MINUSP) {

+      read->querybreak1_strand = '-';

+      read->querybreak1_chrpos = chrpos_low;

+    } else {

+      read->querybreak1_strand = '+';

+      read->querybreak1_chrpos = chrpos_high;

+    }

+  }

+

+  return;

+}

+

+

 static T

 Read_create_bam (bool *completep, T read, Bamline_T bamline) {

   unsigned int flag;

@@ -220,6 +343,7 @@ Read_create_bam (bool *completep, T read, Bamline_T bamline) {

   char strand;

   int type;

   int cigar_querylength, mlength, querypos, i;

+  /* bool insidep; */

   flag = Bamline_flag(bamline);

   shortread = Bamline_read(bamline);

@@ -235,6 +359,9 @@ Read_create_bam (bool *completep, T read, Bamline_T bamline) {

       strcpy(read->quality_string,quality_string);

     }

     read->splicedp = true;

+    read->querybreak1_chrpos = read->querybreak2_chrpos = 0;

+    Read_add_chrinfo_querystart(read,flag,Bamline_chr(bamline),Bamline_chrpos_low(bamline),Bamline_chrpos_high(bamline));

+    Read_add_chrinfo_queryend(read,flag,Bamline_chr(bamline),Bamline_chrpos_low(bamline),Bamline_chrpos_high(bamline));

     *completep = true;

     return read;

   }

@@ -244,6 +371,7 @@ Read_create_bam (bool *completep, T read, Bamline_T bamline) {

   cigar_querylength = Bamline_cigar_querylength(bamline);

   if (read == NULL) {

+    /* Initialize */

     read = (T) MALLOC(sizeof(struct T));

     read->queryseq = (char *) CALLOC(cigar_querylength+1,sizeof(char));

     for (i = 0; i < cigar_querylength; i++) {

@@ -255,24 +383,42 @@ Read_create_bam (bool *completep, T read, Bamline_T bamline) {

       read->quality_string = (char *) CALLOC(cigar_querylength+1,sizeof(char));

     }

     read->splicedp = false;

+    read->querybreak1_chrpos = read->querybreak2_chrpos = 0;

   }

   if (flag & QUERY_MINUSP) {

     strand = '-';

     querypos = cigar_querylength - 1;

+    if (Intlist_head(types) == 'H') {

+      Read_add_chrinfo_querybreak1(read,flag,Bamline_chr(bamline),Bamline_chrpos_low(bamline),Bamline_chrpos_high(bamline));

+    } else {

+      Read_add_chrinfo_queryend(read,flag,Bamline_chr(bamline),Bamline_chrpos_low(bamline),Bamline_chrpos_high(bamline));

+    }

   } else {

     strand = '+';

     querypos = 0;

+    if (Intlist_head(types) == 'H') {

+      Read_add_chrinfo_querybreak2(read,flag,Bamline_chr(bamline),Bamline_chrpos_low(bamline),Bamline_chrpos_high(bamline));

+    } else {

+      Read_add_chrinfo_querystart(read,flag,Bamline_chr(bamline),Bamline_chrpos_low(bamline),Bamline_chrpos_high(bamline));

+    }

   }

+

+

   p = shortread;

   q = quality_string;

+  /* insidep = true; */

   for (a = types, b = npositions; a != NULL; a = Intlist_next(a), b = Uintlist_next(b)) {

     type = Intlist_head(a);

     mlength = Uintlist_head(b);

     if (type == 'H') {

-      querypos += ((strand == '+') ? +mlength : -mlength);

+      if (strand == '+') {

+	querypos += mlength;

+      } else {

+	querypos -= mlength;

+      }

     } else if (type == 'N') {

       /* Ignore */

@@ -328,6 +474,20 @@ Read_create_bam (bool *completep, T read, Bamline_T bamline) {

     }

   }

+  if (type == 'H') {

+    if (flag & QUERY_MINUSP) {

+      Read_add_chrinfo_querybreak2(read,flag,Bamline_chr(bamline),Bamline_chrpos_low(bamline),Bamline_chrpos_high(bamline));

+    } else {

+      Read_add_chrinfo_querybreak1(read,flag,Bamline_chr(bamline),Bamline_chrpos_low(bamline),Bamline_chrpos_high(bamline));

+    }

+  } else {

+    if (flag & QUERY_MINUSP) {

+      Read_add_chrinfo_querystart(read,flag,Bamline_chr(bamline),Bamline_chrpos_low(bamline),Bamline_chrpos_high(bamline));

+    } else {

+      Read_add_chrinfo_queryend(read,flag,Bamline_chr(bamline),Bamline_chrpos_low(bamline),Bamline_chrpos_high(bamline));

+    }

+  }

+

   *completep = true;

   for (i = 0; i < cigar_querylength; i++) {

     if (read->queryseq[i] == '*') {

@@ -335,6 +495,16 @@ Read_create_bam (bool *completep, T read, Bamline_T bamline) {

     }

   }

+#if 0

+  if (insidep == true && Bamline_paired_read_p(bamline) == true) {

+    if (Bamline_lowend_p(bamline) == true) {

+      read->lowp = true;

+    } else {

+      read->highp = true;

+    }

+  }

+#endif

+

   return read;

 }

@@ -349,14 +519,20 @@ Read_create_sam (bool *completep, T read, char *line) {

   int type;

   int cigar_querylength, mlength, querypos, i;

-  Genomicpos_T chrpos_low;

-  char *chr, *acc, *cigar, *auxinfo;

+  Genomicpos_T chrpos_low, chrpos_high, mate_chrpos_low;

+  char *chr, *mate_chr, *acc, *cigar, *auxinfo;

   int mapq;

   int readlength;

+  /* bool insidep; */

-  auxinfo = Samread_parse_line(&acc,&flag,&mapq,&chr,&chrpos_low,&cigar,&readlength,&shortread,&quality_string,line);

+  auxinfo = Samread_parse_line(&acc,&flag,&mapq,&chr,&chrpos_low,&cigar,

+			       &mate_chr,&mate_chrpos_low,&readlength,

+			       &shortread,&quality_string,line);

+  types = Samread_parse_cigar(&npositions,&cigar_querylength,cigar);

+  chrpos_high = Samread_chrpos_high(types,npositions,chrpos_low);

+  FREE(mate_chr);

+  FREE(cigar);

   FREE(acc);

-  FREE(chr);

   if (flag & QUERY_UNMAPPED) {

     read = (T) MALLOC(sizeof(struct T));

@@ -369,19 +545,24 @@ Read_create_sam (bool *completep, T read, char *line) {

       strcpy(read->quality_string,quality_string);

     }

     read->splicedp = true;

+    read->querybreak1_chrpos = read->querybreak2_chrpos = 0;

+    Read_add_chrinfo_querystart(read,flag,chr,chrpos_low,chrpos_high);

+    Read_add_chrinfo_queryend(read,flag,chr,chrpos_low,chrpos_high);

     *completep = true;

+    FREE(chr);

     FREE(cigar);

     FREE(quality_string);

     FREE(shortread);

+    Uintlist_free(&npositions);

+    Intlist_free(&types);

     return read;

   }

-  types = Samread_parse_cigar(&npositions,&cigar_querylength,cigar);

-  FREE(cigar);

   if (read == NULL) {

+    /* Initialize */

     read = (T) MALLOC(sizeof(struct T));

     read->queryseq = (char *) CALLOC(cigar_querylength+1,sizeof(char));

     for (i = 0; i < cigar_querylength; i++) {

@@ -393,24 +574,52 @@ Read_create_sam (bool *completep, T read, char *line) {

       read->quality_string = (char *) CALLOC(cigar_querylength+1,sizeof(char));

     }

     read->splicedp = false;

+    read->querybreak1_chrpos = read->querybreak2_chrpos = 0;

   }

   if (flag & QUERY_MINUSP) {

     strand = '-';

     querypos = cigar_querylength - 1;

+    if (Intlist_head(types) == 'H') {

+      Read_add_chrinfo_querybreak1(read,flag,chr,chrpos_low,chrpos_high);

+    } else {

+      Read_add_chrinfo_queryend(read,flag,chr,chrpos_low,chrpos_high);

+    }

   } else {

     strand = '+';

     querypos = 0;

+    if (Intlist_head(types) == 'H') {

+      Read_add_chrinfo_querybreak2(read,flag,chr,chrpos_low,chrpos_high);