Merge pull request #42 from wtsi-npg/devel

merge devel to master, preparation for release 0.9.4

ChangeLog

@@ -1,3 +1,5 @@
+    Implemented 'bambi chrsplit' for SplitBamByChromosome
+    Implemented 'bambi select' for alignment filtering
 
 	[0.9.3]
 	changes to metrics file output for backwards compatibility (tag#0: no tag name, and perfect_matches pf_perfect matches set to 0)

Makefile.am

@@ -1,5 +1,13 @@
 ACLOCAL_AMFLAGS= -I m4
-TESTS = test/t_bclfile test/t_decode test/t_filterfile test/t_posfile test/t_i2b
+TESTS = test/t_array \
+        test/t_bamit \
+        test/t_chrsplit \
+        test/t_select \
+        test/t_bclfile \
+        test/t_decode \
+        test/t_filterfile \
+        test/t_posfile \
+        test/t_i2b
 SUBDIRS = src test
 dist_doc_DATA = README.md LICENSE
 AM_COLOR_TESTS=always

README.md

@@ -6,15 +6,15 @@ to install:
 
 clone into a local directory (eg 'bambi')
 
-cd bambi
-autoreconf -i
-configure --prefix=<place to install> --with-htslib=<directory containing a htslib release>
+`cd bambi`
 
-(eg configure --prefix=/usr/local/bambi --with-htslib=/usr/local/htslib/1.3)
+`autoreconf -i`
 
-make
+`configure --prefix=<place to install> --with-htslib=<directory containing a htslib release>`
 
-make check (to run tests)
+(eg configure --prefix=/usr/local/bambi --with-htslib=/usr/local/htslib/1.3)
 
+`make`
 
+`make check` (to run tests)
 

configure.ac

@@ -1,6 +1,5 @@
 AC_INIT([bambi], m4_esyscmd_s([git describe --dirty --always --tags]), [[email protected]])
 AC_ARG_VAR(HTSDIR,Directory to look for hts)
-#AM_INIT_AUTOMAKE([-Wall -Werror foreign subdir-objects tar-pax no-dependencies])
 AM_INIT_AUTOMAKE([-Wall -Werror foreign tar-pax no-dependencies])
 AC_PROG_CC
 

src/Makefile.am

@@ -1,9 +1,11 @@
-AM_CPPFLAGS = $(HTSLIB_CPPFLAGS)
+AM_CPPFLAGS = $(HTSLIB_CPPFLAGS) -std=gnu99
 AM_LDFLAGS = -rdynamic $(HTSLIB_LDFLAGS)
 
 bin_PROGRAMS = bambi 
-bambi_SOURCES = bambi.c bclfile.c  decode.c  filterfile.c  hts_addendum.c  i2b.c  posfile.c array.c \
-                bambi.h bclfile.h            filterfile.h  hts_addendum.h         posfile.h array.h
+bambi_SOURCES = bambi.c bclfile.c  decode.c  filterfile.c  hts_addendum.c  i2b.c  posfile.c array.c select.c \
+                bambi.h bclfile.h            filterfile.h  hts_addendum.h         posfile.h array.h \
+                parse.c bamit.c    chrsplit.c\
+                parse.h bamit.h
 nobase_include_HEADERS = cram/cram_samtools.h cram/pooled_alloc.h cram/sam_header.h cram/string_alloc.h
 bambi_CFLAGS = -I/usr/include/libxml2
 bambi_LDADD = $(HTSLIB_HOME)/lib/libhts.a -ldl -lxml2 -lz -llzma -lbz2 -lpthread

src/array.c

@@ -17,6 +17,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
 #include "array.h"
+#include <assert.h>
 
 /*
  * integer array functions
@@ -58,6 +59,34 @@ ia_t *ia_init(int max)
     return ia;
 }
 
+char *ia_join(ia_t *ia, char *delim)
+{
+    int m = 64;
+    char *s = calloc(m,1);
+    char *a = calloc(64,1);
+
+    for (int n=0; n < ia->end; n++) {
+        while (strlen(s)+strlen(delim)>=m) { m *= 2; s = realloc(s,m); }
+        if (n) strcat(s,delim);
+        sprintf(a,"%d",ia->entries[n]);
+        assert(strlen(a)<64);
+        while (strlen(s)+strlen(a)>=m) { m *= 2; s = realloc(s,m); }
+        strcat(s,a);
+        assert(strlen(s) < m);
+    }
+    free(a);
+    return s;
+}
+
+int ia_sum(ia_t *ia)
+{
+    int sum=0;
+    for (int n=0; n < ia->end; n++) {
+        sum += ia->entries[n];
+    }
+    return sum;
+}
+
 
 /*
  * generic arrays
@@ -86,12 +115,29 @@ void va_push(va_t *va, void *ent)
 
 void va_free(va_t *va)
 {
-    int n;
     if (!va) return;
-    for (n=0; n < va->end; n++) {
+    for (int n=0; n < va->end; n++) {
         va->free_entry(va->entries[n]);
     }
     free(va->entries);
     free(va);
 }
 
+char *va_join(va_t *va, char *delim)
+{
+    int m = 64;
+    char *s = calloc(m,1);
+    char *a;
+
+    for (int n=0; n < va->end; n++) {
+        while (strlen(s)+strlen(delim)>=m) { m *= 2; s = realloc(s,m); }
+        if (n) strcat(s,delim);
+        a = va->entries[n];
+        while (strlen(s)+strlen(a)>=m) { m *= 2; s = realloc(s,m); }
+        strcat(s,a);
+        assert(strlen(s) < m);
+    }
+    return s;
+}
+
+

src/array.h

@@ -21,6 +21,8 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 #include <stdbool.h>
 #include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
 
 /*
  * integer array functions
@@ -37,6 +39,8 @@ int ia_compare(const void *ia1, const void *ia2);
 void ia_sort(ia_t *ia);
 void ia_push(ia_t *ia, int i);
 void ia_free(ia_t *ia);
+char *ia_join(ia_t *ia, char *delim);
+int ia_sum(ia_t *ia);
 static inline bool ia_isEmpty(ia_t *ia) { return (ia->end == 0); }
 
 
@@ -55,6 +59,7 @@ va_t *va_init(int max, void(*free_entry)(void*));
 void va_push(va_t *va, void *ent);
 void va_free(va_t *va);
 static inline bool va_isEmpty(va_t *va) { return va->end == 0; }
+char *va_join(va_t *va, char *delim);
 
 #endif
 

src/bambi.c

@@ -31,7 +31,8 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 int main_decode(int argc, char *argv[]);
 int main_i2b(int argc, char *argv[]);
-int main_afilter(int argc, char *argv[]);
+int main_select(int argc, char *argv[]);
+int main_chrsplit(int argc, char *argv[]);
 
 const char *bambi_version()
 {
@@ -80,7 +81,10 @@ static void usage(FILE *fp)
 "Commands:\n"
 "     decode         decode a multiplexed SAM/BAM/CRAM file by read groups\n"
 "     i2b            converts illumina files to SAM/BAM/CRAM files\n"
-"     afilter        alignment filter\n"
+"     select         select reads by alignment\n"
+"     chrsplit       split reads by chromosome\n"
+"\n"
+"bambi <command> for help on a particular command\n"
 "\n");
 }
 
@@ -101,6 +105,8 @@ int main(int argc, char *argv[])
     int ret = 0;
          if (strcmp(argv[1], "decode") == 0)    ret = main_decode(argc-1, argv+1);
     else if (strcmp(argv[1], "i2b") == 0)       ret = main_i2b(argc-1, argv+1);
+    else if (strcmp(argv[1], "select") == 0)    ret = main_select(argc-1, argv+1);
+    else if (strcmp(argv[1], "chrsplit") == 0)  ret = main_chrsplit(argc-1, argv+1);
     else if (strcmp(argv[1], "--version") == 0) {
         printf( "bambi %s\n"
                 "Using htslib %s\n"

src/bambi.h

@@ -21,9 +21,14 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #ifndef __BAMBI_H__
 #define __BAMBI_H__
 
+#define INDEX_SEPARATOR "-"
+
 #include "config.h"
 #include "hts_addendum.h"
+#include "array.h"
+
 const char *bambi_version(void);
+void parse_tags(va_t *tags, char *arg);
 
 #endif
 

src/bamit.c

@@ -0,0 +1,65 @@
+/*  bamit.c -- iterator for reading BAM files one record at a time.
+
+    Copyright (C) 2016 Genome Research Ltd.
+
+    Author: Jennifer Liddle <[email protected]>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Affero General Public License as published
+by the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "bamit.h"
+
+BAMit_t *BAMit_init(samFile *f, bam_hdr_t *h)
+{
+    int r;
+    BAMit_t *bit = calloc(1,sizeof(BAMit_t));
+    bit->f = f;
+    bit->h = h;
+    bit->rec = bam_init1();
+    bit->nextRec = bam_init1();
+    if (f->is_write == 0) r = sam_read1(bit->f, bit->h, bit->nextRec);
+    return bit;
+}
+
+void BAMit_free(void *ptr)
+{
+    BAMit_t *bit = (BAMit_t *)ptr;
+    if (!bit) return;
+    if (bit->f) hts_close(bit->f);
+    if (bit->h) bam_hdr_destroy(bit->h);
+    if (bit->rec) bam_destroy1(bit->rec);
+    if (bit->nextRec) bam_destroy1(bit->nextRec);
+    free(bit);
+}
+
+bam1_t *BAMit_next(BAMit_t *bit)
+{
+    if (!bit->nextRec) return NULL;
+    bam_copy1(bit->rec,bit->nextRec);
+    int r = sam_read1(bit->f, bit->h, bit->nextRec);
+    if (r<0) { bam_destroy1(bit->nextRec); bit->nextRec = NULL; }
+    return bit->rec;
+}
+
+bam1_t *BAMit_peek(BAMit_t *bit)
+{
+    return bit->nextRec;
+}
+
+bool BAMit_hasnext(BAMit_t *bit)
+{
+    return (bit->nextRec != NULL);
+}
+
+
+
+

src/bamit.h

@@ -0,0 +1,44 @@
+/*  bamit.h -- BAM Iterator - for reading BAM files one record at a time.
+
+    Copyright (C) 2016 Genome Research Ltd.
+
+    Author: Jennifer Liddle <[email protected]>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Affero General Public License as published
+by the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __BAMIT_H__
+#define __BAMIT_H__
+
+#include <stdbool.h>
+#include <stdlib.h>
+#include "htslib/sam.h"
+
+/*
+ * iterator structure
+ */
+
+typedef struct {
+    samFile *f;
+    bam_hdr_t *h;
+    bam1_t *rec;
+    bam1_t *nextRec;
+} BAMit_t;
+
+BAMit_t *BAMit_init(samFile *f, bam_hdr_t *h);
+void BAMit_free(void *bit);
+bam1_t *BAMit_next(BAMit_t *bit);
+bam1_t *BAMit_peek(BAMit_t *bit);
+bool BAMit_hasnext(BAMit_t *bit);
+
+#endif
+