* put version number in bam.h

 * write version to BCF
 * in phase, change the default -q to 37
 * output a little more information during phasing

diff --git a/Makefile b/Makefile
index 80f7a31b94f5d2a30909a0cac65a3db8c0f8462f..d557520dd5b0c60b7ab4763e34939d3031b2acb1 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
 CC=                    gcc
 CFLAGS=                -g -Wall -O2 #-m64 #-arch ppc
-DFLAGS=                -D_FILE_OFFSET_BITS=64 -D_USE_KNETFILE -D_CURSES_LIB=1
+DFLAGS=                -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE -D_USE_KNETFILE -D_CURSES_LIB=1
 KNETFILE_O=    knetfile.o
 LOBJS=         bgzf.o kstring.o bam_aux.o bam.o bam_import.o sam.o bam_index.o \
                        bam_pileup.o bam_lpileup.o bam_md.o glf.o razf.o faidx.o \

diff --git a/bam.h b/bam.h
index eef2ea9853f25966c5a40c45a86272b78a507fb2..752d0bda7dea595c3c9014968fee9cd02057662b 100644 (file)
--- a/bam.h
+++ b/bam.h
@@ -40,6 +40,8 @@
   @copyright Genome Research Ltd.
  */
 
+#define BAM_VERSION "0.1.12-r921:124"
+
 #include <stdint.h>
 #include <stdlib.h>
 #include <string.h>

diff --git a/bam_plcmd.c b/bam_plcmd.c
index 2c333e8c668cb1a48464637d2e56e2e3984d7269..bba62cbcdcabe52cbdba71ab71a6b2e5681aa032 100644 (file)
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -701,7 +701,8 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                bh->l_smpl = s.l;
                bh->sname = malloc(s.l);
                memcpy(bh->sname, s.s, s.l);
-               bh->l_txt = 0;
+               bh->txt = malloc(strlen(BAM_VERSION) + 64);
+               bh->l_txt = 1 + sprintf(bh->txt, "##samtoolsVersion=%s\n", BAM_VERSION);
                free(s.s);
                bcf_hdr_sync(bh);
                bcf_hdr_write(bp, bh);

diff --git a/bamtk.c b/bamtk.c
index 03e413dc966f6f582014832db8aca307e5725761..1d1151973fb08efbbb6d642480712fbc6b3b1563 100644 (file)
--- a/bamtk.c
+++ b/bamtk.c
@@ -8,10 +8,6 @@
 #include "knetfile.h"
 #endif
 
-#ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.1.12-11 (r921:118)"
-#endif
-
 int bam_taf2baf(int argc, char *argv[]);
 int bam_pileup(int argc, char *argv[]);
 int bam_mpileup(int argc, char *argv[]);
@@ -77,7 +73,7 @@ static int usage()
 {
        fprintf(stderr, "\n");
        fprintf(stderr, "Program: samtools (Tools for alignments in the SAM format)\n");
-       fprintf(stderr, "Version: %s\n\n", PACKAGE_VERSION);
+       fprintf(stderr, "Version: %s\n\n", BAM_VERSION);
        fprintf(stderr, "Usage:   samtools <command> [options]\n\n");
        fprintf(stderr, "Command: view        SAM<->BAM conversion\n");
        fprintf(stderr, "         sort        sort alignment file\n");

diff --git a/phase.c b/phase.c
index 7f069584b9fe4d1cee05514ae7286a08cfd8a2db..5fc4ce77d94ed3046082abf7a9525caf3362f7c3 100644 (file)
--- a/phase.c
+++ b/phase.c
@@ -35,6 +35,7 @@ typedef struct {
        int8_t seq[MAX_VARS]; // TODO: change to dynamic memory allocation!
        int vpos, beg, end;
        uint32_t vlen:16, single:1, flip:1, phase:1, phased:1;
+       uint32_t in:16, out:16; // in-phase and out-phase
 } frag_t, *frag_p;
 
 #define rseq_lt(a,b) ((a)->vpos < (b)->vpos)
@@ -174,7 +175,8 @@ static uint64_t *fragphase(int vpos, const int8_t *path, nseq_t *hash, int flip)
                                ++c[f->seq[i] == path[f->vpos + i] + 1? 0 : 1];
                        }
                        f->phase = c[0] > c[1]? 0 : 1;
-                       f->phased = c[0] == c[1]? 0 : 1;
+                       f->in = c[f->phase]; f->out = c[1 - f->phase];
+                       if (f->in && f->out && f->in <= f->out + 1) f->phased = 0;
                        // fix chimera
                        f->flip = 0;
                        if (flip && c[0] >= 3 && c[1] >= 3) {
@@ -320,7 +322,7 @@ static void dump_aln(phaseg_t *g, int min_pos, const nseq_t *hash)
                else {
                        frag_t *f = &kh_val(hash, k);
                        if (f->phased && f->flip) which = 2;
-                       else if (f->phased == 0) which = 3;
+                       else if (f->phased == 0) which = 2;
                        else { // phased and not flipped
                                char c = 'Y';
                                which = f->phase;
@@ -395,7 +397,7 @@ static int phase(phaseg_t *g, const char *chr, int vpos, uint64_t *cns, nseq_t *
                int8_t c[2];
                c[0] = (cns[i]&0xffff)>>2 == 0? 4 : (cns[i]&3);
                c[1] = (cns[i]>>16&0xffff)>>2 == 0? 4 : (cns[i]>>16&3);
-               printf("M%d\t%s\t%d\t%d\t%c\t%c\t%d\t%d\t%d\t%d\t%d\n", sitemask[i]+1, chr, (int)(cns[0]>>32), (int)(cns[i]>>32) + 1, "ACGTX"[c[path[i]]], "ACGTX"[c[1-path[i]]],
+               printf("M%d\t%s\t%d\t%d\t%c\t%c\t%d\t%d\t%d\t%d\t%d\n", sitemask[i]+1, chr, (int)(cns[0]>>32) + 1, (int)(cns[i]>>32) + 1, "ACGTX"[c[path[i]]], "ACGTX"[c[1-path[i]]],
                        i + g->vpos_shift + 1, (int)(x&0xffff), (int)(x>>16&0xffff), (int)(x>>32&0xffff), (int)(x>>48&0xffff));
        }
        free(path); free(pcnt); free(regmask); free(sitemask);
@@ -413,7 +415,7 @@ static int phase(phaseg_t *g, const char *chr, int vpos, uint64_t *cns, nseq_t *
                        if (f->seq[j] == 0) putchar('N');
                        else putchar("ACGT"[f->seq[j] == 1? (c&3) : (c>>16&3)]);
                }
-               printf("\t*\tYP:i:%d\tYF:i:%d\n", f->phase, f->flip);
+               printf("\t*\tYP:i:%d\tYF:i:%d\tYI:i:%d\tYO:i:%d\tYS:i:%d\n", f->phase, f->flip, f->in, f->out, f->beg+1);
        }
        free(seqs);
        printf("//\n");
@@ -445,6 +447,7 @@ static int readaln(void *data, bam1_t *b)
        phaseg_t *g = (phaseg_t*)data;
        int ret;
        ret = bam_read1(g->fp, b);
+       if (ret < 0) return ret;
        if (!(b->core.flag & (BAM_FUNMAP|BAM_FSECONDARY|BAM_FQCFAIL|BAM_FDUP)) && g->pre) {
                if (g->n == g->m) {
                        g->m = g->m? g->m<<1 : 16;
@@ -515,7 +518,7 @@ int main_phase(int argc, char *argv[])
 
        memset(&g, 0, sizeof(phaseg_t));
        g.flag = FLAG_FIX_CHIMERA;
-       g.min_varLOD = 40; g.k = 13; g.min_baseQ = 13; g.max_depth = 256;
+       g.min_varLOD = 37; g.k = 13; g.min_baseQ = 13; g.max_depth = 256;
        while ((c = getopt(argc, argv, "Q:eFq:k:b:l:D:")) >= 0) {
                switch (c) {
                        case 'D': g.max_depth = atoi(optarg); break;