* samtools-0.1.9-8 (r807)

 * collect indel candidates only from specified platforms (@RG-PL)
 * merge varFilter and filter4vcf in vcfutils.pl

diff --git a/Makefile b/Makefile
index eeae213875b6b084592d7f0858452915f2a53695..965491a9fa518bc990328b34699863c418554cb3 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -64,6 +64,7 @@ glf.o:glf.h
 sam_header.o:sam_header.h khash.h
 bcf.o:bcftools/bcf.h
 bam2bcf.o:bam2bcf.h errmod.h bcftools/bcf.h
 sam_header.o:sam_header.h khash.h
 bcf.o:bcftools/bcf.h
 bam2bcf.o:bam2bcf.h errmod.h bcftools/bcf.h

+bam2bcf_indel.o:bam2bcf.h

 errmod.o:errmod.h
 
 faidx.o:faidx.h razf.h khash.h
 errmod.o:errmod.h
 
 faidx.o:faidx.h razf.h khash.h

diff --git a/bam2bcf.h b/bam2bcf.h
index 570990e4174caec6220c130e0c343dec1d6dd562..e0673e5e868b3ab4babcd0726f76ff301c42e9be 100644 (file)
--- a/bam2bcf.h
+++ b/bam2bcf.h
@@ -17,6 +17,7 @@ typedef struct __bcf_callaux_t {
        char *inscns;
        uint16_t *bases;
        errmod_t *e;
        char *inscns;
        uint16_t *bases;
        errmod_t *e;

+       void *rghash;

 } bcf_callaux_t;
 
 typedef struct {
 } bcf_callaux_t;
 
 typedef struct {


@@ -42,7 +43,8 @@ extern "C" {
        int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, bcf_call_t *call);
        int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bcr, int is_SP,
                                         const bcf_callaux_t *bca, const char *ref);
        int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, bcf_call_t *call);
        int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bcr, int is_SP,
                                         const bcf_callaux_t *bca, const char *ref);

-       int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref);
+       int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref,
+                                                 const void *rghash);

 
 #ifdef __cplusplus
 }
 
 #ifdef __cplusplus
 }

diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index de9df3e6c2fb9042b4287ac60e828224aaa651a2..2d0b5c6f8a66a4c5c006dbec560718ab1f7a0377 100644 (file)
--- a/bam2bcf_indel.c
+++ b/bam2bcf_indel.c
@@ -1,13 +1,61 @@
 #include <assert.h>
 #include <ctype.h>
 #include <assert.h>
 #include <ctype.h>

+#include <string.h>

 #include "bam.h"
 #include "bam2bcf.h"
 #include "ksort.h"
 #include "kaln.h"
 #include "bam.h"
 #include "bam2bcf.h"
 #include "ksort.h"
 #include "kaln.h"

+#include "khash.h"
+KHASH_SET_INIT_STR(rg)

 
 #define MINUS_CONST 0x10000000
 #define INDEL_WINDOW_SIZE 50
 
 
 #define MINUS_CONST 0x10000000
 #define INDEL_WINDOW_SIZE 50
 

+void *bcf_call_add_rg(void *_hash, const char *hdtext, const char *list)
+{
+       const char *s, *p, *q, *r, *t;
+       khash_t(rg) *hash;
+       if (list == 0 || hdtext == 0) return _hash;
+       if (_hash == 0) _hash = kh_init(rg);
+       hash = (khash_t(rg)*)_hash;
+       if ((s = strstr(hdtext, "@RG\t")) == 0) return hash;
+       do {
+               t = strstr(s + 4, "@RG\t"); // the next @RG
+               if ((p = strstr(s, "\tID:")) != 0) p += 4;
+               if ((q = strstr(s, "\tPL:")) != 0) q += 4;
+               if (p && q && (t == 0 || (p < t && q < t))) { // ID and PL are both present
+                       int lp, lq;
+                       char *x;
+                       for (r = p; *r && *r != '\t' && *r != '\n'; ++r); lp = r - p;
+                       for (r = q; *r && *r != '\t' && *r != '\n'; ++r); lq = r - q;
+                       x = calloc((lp > lq? lp : lq) + 1, 1);
+                       for (r = q; *r && *r != '\t' && *r != '\n'; ++r) x[r-q] = *r;
+                       if (strstr(list, x)) { // insert ID to the hash table
+                               khint_t k;
+                               int ret;
+                               for (r = p; *r && *r != '\t' && *r != '\n'; ++r) x[r-p] = *r;
+                               x[r-p] = 0;
+                               k = kh_get(rg, hash, x);
+                               if (k == kh_end(hash)) k = kh_put(rg, hash, x, &ret);
+                               else free(x);
+                       } else free(x);
+               }
+               s = t;
+       } while (s);
+       return hash;
+}
+
+void bcf_call_del_rghash(void *_hash)
+{
+       khint_t k;
+       khash_t(rg) *hash = (khash_t(rg)*)_hash;
+       if (hash == 0) return;
+       for (k = kh_begin(hash); k < kh_end(hash); ++k)
+               if (kh_exist(hash, k))
+                       free((char*)kh_key(hash, k));
+       kh_destroy(rg, hash);
+}
+

 static int tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos)
 {
        int k, x = c->pos, y = 0, last_y = 0;
 static int tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos)
 {
        int k, x = c->pos, y = 0, last_y = 0;


@@ -58,12 +106,30 @@ static inline int est_indelreg(int pos, const char *ref, int l, char *ins4)
        return max_i - pos;
 }
 
        return max_i - pos;
 }
 

-int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref)
+int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref,
+                                         const void *rghash)

 {
        extern void ks_introsort_uint32_t(int, uint32_t*);
        int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score, N, K, l_run, ref_type;
        char *inscns = 0, *ref2, *query;
 {
        extern void ks_introsort_uint32_t(int, uint32_t*);
        int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score, N, K, l_run, ref_type;
        char *inscns = 0, *ref2, *query;

+       khash_t(rg) *hash = (khash_t(rg)*)rghash;

        if (ref == 0 || bca == 0) return -1;
        if (ref == 0 || bca == 0) return -1;

+       // mark filtered reads
+       if (rghash) {
+               N = 0;
+               for (s = N = 0; s < n; ++s) {
+                       for (i = 0; i < n_plp[s]; ++i) {
+                               bam_pileup1_t *p = plp[s] + i;
+                               const uint8_t *rg = bam_aux_get(p->b, "RG");
+                               p->aux = 1; // filtered by default
+                               if (rg) {
+                                       khint_t k = kh_get(rg, hash, (const char*)(rg + 1));
+                                       if (k != kh_end(hash)) p->aux = 0, ++N; // not filtered
+                               }
+                       }
+               }
+               if (N == 0) return -1; // no reads left
+       }

        // determine if there is a gap
        for (s = N = 0; s < n; ++s) {
                for (i = 0; i < n_plp[s]; ++i)
        // determine if there is a gap
        for (s = N = 0; s < n; ++s) {
                for (i = 0; i < n_plp[s]; ++i)


@@ -81,7 +147,7 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                for (s = 0; s < n; ++s) {
                        for (i = 0; i < n_plp[s]; ++i) {
                                const bam_pileup1_t *p = plp[s] + i;
                for (s = 0; s < n; ++s) {
                        for (i = 0; i < n_plp[s]; ++i) {
                                const bam_pileup1_t *p = plp[s] + i;

-                               if (p->indel != 0)
+                               if (p->indel != 0 && (rghash == 0 || p->aux == 0))

                                        aux[m++] = MINUS_CONST + p->indel;
                                j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b));
                                if (j > max_rd_len) max_rd_len = j;
                                        aux[m++] = MINUS_CONST + p->indel;
                                j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b));
                                if (j > max_rd_len) max_rd_len = j;


@@ -91,7 +157,9 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                // squeeze out identical types
                for (i = 1, n_types = 1; i < m; ++i)
                        if (aux[i] != aux[i-1]) ++n_types;
                // squeeze out identical types
                for (i = 1, n_types = 1; i < m; ++i)
                        if (aux[i] != aux[i-1]) ++n_types;

-               assert(n_types > 1); // there must at least one type of non-reference indel
+               if (n_types == 1) { // no indels
+                       free(aux); return -1;
+               }

                types = (int*)calloc(n_types, sizeof(int));
                t = 0;
                types[t++] = aux[0] - MINUS_CONST; 
                types = (int*)calloc(n_types, sizeof(int));
                t = 0;
                types[t++] = aux[0] - MINUS_CONST; 

diff --git a/bam_plcmd.c b/bam_plcmd.c
index 54b8f91722b2f16ade955231e075764b3161d2d5..0fbc9303c25744ebde8d1b0e269c539d5c4f5325 100644 (file)
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -535,7 +535,7 @@ int bam_pileup(int argc, char *argv[])
 
 typedef struct {
        int max_mq, min_mq, flag, min_baseQ, capQ_thres, max_depth;
 
 typedef struct {
        int max_mq, min_mq, flag, min_baseQ, capQ_thres, max_depth;

-       char *reg, *fn_pos;
+       char *reg, *fn_pos, *pl_list;

        faidx_t *fai;
        kh_64_t *hash;
 } mplp_conf_t;
        faidx_t *fai;
        kh_64_t *hash;
 } mplp_conf_t;


@@ -602,6 +602,8 @@ static void group_smpl(mplp_pileup_t *m, bam_sample_t *sm, kstring_t *buf,
 
 static int mpileup(mplp_conf_t *conf, int n, char **fn)
 {
 
 static int mpileup(mplp_conf_t *conf, int n, char **fn)
 {

+       extern void *bcf_call_add_rg(void *rghash, const char *hdtext, const char *list);
+       extern void bcf_call_del_rghash(void *rghash);

        mplp_aux_t **data;
        int i, tid, pos, *n_plp, beg0 = 0, end0 = 1u<<29, ref_len, ref_tid, max_depth;
        const bam_pileup1_t **plp;
        mplp_aux_t **data;
        int i, tid, pos, *n_plp, beg0 = 0, end0 = 1u<<29, ref_len, ref_tid, max_depth;
        const bam_pileup1_t **plp;


@@ -609,6 +611,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
        bam_header_t *h = 0;
        char *ref;
        khash_t(64) *hash = 0;
        bam_header_t *h = 0;
        char *ref;
        khash_t(64) *hash = 0;

+       void *rghash = 0;

 
        bcf_callaux_t *bca = 0;
        bcf_callret1_t *bcr = 0;
 
        bcf_callaux_t *bca = 0;
        bcf_callret1_t *bcr = 0;


@@ -638,6 +641,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                data[i]->fp = strcmp(fn[i], "-") == 0? bam_dopen(fileno(stdin), "r") : bam_open(fn[i], "r");
                h_tmp = bam_header_read(data[i]->fp);
                bam_smpl_add(sm, fn[i], h_tmp->text);
                data[i]->fp = strcmp(fn[i], "-") == 0? bam_dopen(fileno(stdin), "r") : bam_open(fn[i], "r");
                h_tmp = bam_header_read(data[i]->fp);
                bam_smpl_add(sm, fn[i], h_tmp->text);

+               rghash = bcf_call_add_rg(rghash, h_tmp->text, conf->pl_list);

                if (conf->reg) {
                        int beg, end;
                        bam_index_t *idx;
                if (conf->reg) {
                        int beg, end;
                        bam_index_t *idx;


@@ -693,6 +697,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                bcf_hdr_write(bp, bh);
                bca = bcf_call_init(-1., conf->min_baseQ);
                bcr = calloc(sm->n, sizeof(bcf_callret1_t));
                bcf_hdr_write(bp, bh);
                bca = bcf_call_init(-1., conf->min_baseQ);
                bcr = calloc(sm->n, sizeof(bcf_callret1_t));

+               bca->rghash = rghash;

        }
        ref_tid = -1; ref = 0;
        iter = bam_mplp_init(n, mplp_func, (void**)data);
        }
        ref_tid = -1; ref = 0;
        iter = bam_mplp_init(n, mplp_func, (void**)data);


@@ -731,7 +736,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                        bcf_write(bp, bh, b);
                        bcf_destroy(b);
                        // call indels
                        bcf_write(bp, bh, b);
                        bcf_destroy(b);
                        // call indels

-                       if (bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref) >= 0) {
+                       if (bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref, rghash) >= 0) {

                                for (i = 0; i < gplp.n; ++i)
                                        bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], -1, bca, bcr + i);
                                bcf_call_combine(gplp.n, bcr, -1, &bc);
                                for (i = 0; i < gplp.n; ++i)
                                        bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], -1, bca, bcr + i);
                                bcf_call_combine(gplp.n, bcr, -1, &bc);


@@ -767,6 +772,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
        bam_smpl_destroy(sm); free(buf.s);
        for (i = 0; i < gplp.n; ++i) free(gplp.plp[i]);
        free(gplp.plp); free(gplp.n_plp); free(gplp.m_plp);
        bam_smpl_destroy(sm); free(buf.s);
        for (i = 0; i < gplp.n; ++i) free(gplp.plp[i]);
        free(gplp.plp); free(gplp.n_plp); free(gplp.m_plp);

+       bcf_call_del_rghash(rghash);

        if (hash) { // free the hash table
                khint_t k;
                for (k = kh_begin(hash); k < kh_end(hash); ++k)
        if (hash) { // free the hash table
                khint_t k;
                for (k = kh_begin(hash); k < kh_end(hash); ++k)


@@ -847,7 +853,7 @@ int bam_mpileup(int argc, char *argv[])
        mplp.capQ_thres = 0;
        mplp.max_depth = 250;
        mplp.flag = MPLP_NO_ORPHAN | MPLP_REALN;
        mplp.capQ_thres = 0;
        mplp.max_depth = 250;
        mplp.flag = MPLP_NO_ORPHAN | MPLP_REALN;

-       while ((c = getopt(argc, argv, "gf:r:l:M:q:Q:uaORC:BDSd:b:")) >= 0) {
+       while ((c = getopt(argc, argv, "gf:r:l:M:q:Q:uaORC:BDSd:b:P:")) >= 0) {

                switch (c) {
                case 'f':
                        mplp.fai = fai_load(optarg);
                switch (c) {
                case 'f':
                        mplp.fai = fai_load(optarg);


@@ -856,6 +862,7 @@ int bam_mpileup(int argc, char *argv[])
                case 'd': mplp.max_depth = atoi(optarg); break;
                case 'r': mplp.reg = strdup(optarg); break;
                case 'l': mplp.fn_pos = strdup(optarg); break;
                case 'd': mplp.max_depth = atoi(optarg); break;
                case 'r': mplp.reg = strdup(optarg); break;
                case 'l': mplp.fn_pos = strdup(optarg); break;

+               case 'P': mplp.pl_list = strdup(optarg); break;

                case 'g': mplp.flag |= MPLP_GLF; break;
                case 'u': mplp.flag |= MPLP_NO_COMP | MPLP_GLF; break;
                case 'a': mplp.flag |= MPLP_NO_ORPHAN | MPLP_REALN; break;
                case 'g': mplp.flag |= MPLP_GLF; break;
                case 'u': mplp.flag |= MPLP_NO_COMP | MPLP_GLF; break;
                case 'a': mplp.flag |= MPLP_NO_ORPHAN | MPLP_REALN; break;


@@ -882,6 +889,7 @@ int bam_mpileup(int argc, char *argv[])
                fprintf(stderr, "         -Q INT      min base quality [%d]\n", mplp.min_baseQ);
                fprintf(stderr, "         -q INT      filter out alignment with MQ smaller than INT [%d]\n", mplp.min_mq);
                fprintf(stderr, "         -d INT      max per-sample depth [%d]\n", mplp.max_depth);
                fprintf(stderr, "         -Q INT      min base quality [%d]\n", mplp.min_baseQ);
                fprintf(stderr, "         -q INT      filter out alignment with MQ smaller than INT [%d]\n", mplp.min_mq);
                fprintf(stderr, "         -d INT      max per-sample depth [%d]\n", mplp.max_depth);

+               fprintf(stderr, "         -P STR      comma separated list of platforms for indels [all]\n");

                fprintf(stderr, "         -g          generate BCF output\n");
                fprintf(stderr, "         -u          do not compress BCF output\n");
                fprintf(stderr, "         -B          disable BAQ computation\n");
                fprintf(stderr, "         -g          generate BCF output\n");
                fprintf(stderr, "         -u          do not compress BCF output\n");
                fprintf(stderr, "         -B          disable BAQ computation\n");


@@ -900,7 +908,7 @@ int bam_mpileup(int argc, char *argv[])
     }
     else
            mpileup(&mplp, argc - optind, argv + optind);
     }
     else
            mpileup(&mplp, argc - optind, argv + optind);

-       free(mplp.reg);
+       free(mplp.reg); free(mplp.pl_list);

        if (mplp.fai) fai_destroy(mplp.fai);
        return 0;
 }
        if (mplp.fai) fai_destroy(mplp.fai);
        return 0;
 }

diff --git a/bamtk.c b/bamtk.c
index e76a9d4499acb9404910732ec0ff624ed11c7b1d..d2817b9902a5bccf985bb8360775029bef946267 100644 (file)
--- a/bamtk.c
+++ b/bamtk.c
@@ -9,7 +9,7 @@
 #endif
 
 #ifndef PACKAGE_VERSION
 #endif
 
 #ifndef PACKAGE_VERSION

-#define PACKAGE_VERSION "0.1.9-7 (r804)"
+#define PACKAGE_VERSION "0.1.9-8 (r807)"

 #endif
 
 int bam_taf2baf(int argc, char *argv[]);
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);

diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index 0cf73a30d640f6e0747e2cc4b9af8e4bbfb815a7..c1177eacb9603cc2cb9f93ab27d54b5add1e0779 100755 (executable)
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -10,11 +10,10 @@ use Getopt::Std;
 exit;
 
 sub main {
 exit;
 
 sub main {

-  my $version = '0.1.0';

   &usage if (@ARGV < 1);
   my $command = shift(@ARGV);
   my %func = (subsam=>\&subsam, listsam=>\&listsam, fillac=>\&fillac, qstats=>\&qstats, varFilter=>\&varFilter,
   &usage if (@ARGV < 1);
   my $command = shift(@ARGV);
   my %func = (subsam=>\&subsam, listsam=>\&listsam, fillac=>\&fillac, qstats=>\&qstats, varFilter=>\&varFilter,

-                         hapmap2vcf=>\&hapmap2vcf, ucscsnp2vcf=>\&ucscsnp2vcf, filter4vcf=>\&filter4vcf, ldstats=>\&ldstats);
+                         hapmap2vcf=>\&hapmap2vcf, ucscsnp2vcf=>\&ucscsnp2vcf, filter4vcf=>\&varFilter, ldstats=>\&ldstats);

   die("Unknown command \"$command\".\n") if (!defined($func{$command}));
   &{$func{$command}};
 }
   die("Unknown command \"$command\".\n") if (!defined($func{$command}));
   &{$func{$command}};
 }


@@ -199,8 +198,8 @@ Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions #
 }
 
 sub varFilter {
 }
 
 sub varFilter {

-  my %opts = (d=>1, D=>10000, a=>2, l=>10, Q=>10, s=>100, w=>10, p=>undef, F=>.001);
-  getopts('pd:D:l:Q:w:G:F:a:', \%opts);
+  my %opts = (d=>2, D=>10000, a=>2, W=>10, Q=>10, w=>10, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4);
+  getopts('pd:D:W:Q:w:a:1:2:3:4:', \%opts);

   die(qq/
 Usage:   vcfutils.pl varFilter [options] <in.vcf>
 
   die(qq/
 Usage:   vcfutils.pl varFilter [options] <in.vcf>
 


@@ -209,19 +208,28 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
          -D INT    maximum read depth [$opts{D}]
          -a INT    minimum number of alternate bases [$opts{a}]
          -w INT    SNP within INT bp around a gap to be filtered [$opts{w}]
          -D INT    maximum read depth [$opts{D}]
          -a INT    minimum number of alternate bases [$opts{a}]
          -w INT    SNP within INT bp around a gap to be filtered [$opts{w}]

-         -l INT    window size for filtering adjacent gaps [$opts{l}]
+         -W INT    window size for filtering adjacent gaps [$opts{W}]
+         -1 FLOAT  min P-value for strand bias (given PV4) [$opts{1}]
+         -2 FLOAT  min P-value for baseQ bias [$opts{2}]
+         -3 FLOAT  min P-value for mapQ bias [$opts{3}]
+         -4 FLOAT  min P-value for end distance bias [$opts{4}]

          -p        print filtered variants
          -p        print filtered variants

+
+Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.

 \n/) if (@ARGV == 0 && -t STDIN);
 
   # calculate the window size
 \n/) if (@ARGV == 0 && -t STDIN);
 
   # calculate the window size

-  my ($ol, $ow) = ($opts{l}, $opts{w});
+  my ($ol, $ow) = ($opts{W}, $opts{w});

   my $max_dist = $ol > $ow? $ol : $ow;
   # the core loop
   my @staging; # (indel_filtering_score, flt_tag, indel_span; chr, pos, ...)
   while (<>) {
        my @t = split;
   my $max_dist = $ol > $ow? $ol : $ow;
   # the core loop
   my @staging; # (indel_filtering_score, flt_tag, indel_span; chr, pos, ...)
   while (<>) {
        my @t = split;

-       next if (/^#/);
+    if (/^#/) {
+         print; next;
+       }

        next if ($t[4] eq '.'); # skip non-var sites
        next if ($t[4] eq '.'); # skip non-var sites

+       # check if the site is a SNP

        my $is_snp = 1;
        if (length($t[3]) > 1) {
          $is_snp = 0;
        my $is_snp = 1;
        if (length($t[3]) > 1) {
          $is_snp = 0;


@@ -238,7 +246,6 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
          varFilter_aux(shift(@staging), $opts{p}); # calling a function is a bit slower, not much
        }
        my $flt = 0;
          varFilter_aux(shift(@staging), $opts{p}); # calling a function is a bit slower, not much
        }
        my $flt = 0;

-

        # parse annotations
        my ($dp, $mq, $dp_alt) = (-1, -1, -1);
        if ($t[7] =~ /DP=(\d+)/i) {
        # parse annotations
        my ($dp, $mq, $dp_alt) = (-1, -1, -1);
        if ($t[7] =~ /DP=(\d+)/i) {


@@ -258,6 +265,8 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
        }
        $flt = 4 if ($dp_alt >= 0 && $dp_alt < $opts{a});
        $flt = 1 if ($flt == 0 && $mq >= 0 && $mq < $opts{Q});
        }
        $flt = 4 if ($dp_alt >= 0 && $dp_alt < $opts{a});
        $flt = 1 if ($flt == 0 && $mq >= 0 && $mq < $opts{Q});

+       $flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/
+                                && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));

 
        # site dependent filters
        my ($rlen, $indel_score) = (0, -1); # $indel_score<0 for SNPs
 
        # site dependent filters
        my ($rlen, $indel_score) = (0, -1); # $indel_score<0 for SNPs


@@ -300,46 +309,7 @@ sub varFilter_aux {
   if ($first->[1] == 0) {
        print join("\t", @$first[3 .. @$first-1]), "\n";
   } elsif ($is_print) {
   if ($first->[1] == 0) {
        print join("\t", @$first[3 .. @$first-1]), "\n";
   } elsif ($is_print) {

-       print STDERR join("\t", substr("UQdDaGgsiX", $first->[1], 1), @$first[3 .. @$first-1]), "\n";
-  }
-}
-
-sub filter4vcf {
-  my %opts = (d=>3, D=>2000, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, Q=>10, q=>3);
-  getopts('d:D:1:2:3:4:Q:q:', \%opts);
-  die(qq/
-Usage:   vcfutils.pl filter4vcf [options] <in.vcf>
-
-Options: -d INT     min total depth (given DP or DP4) [$opts{d}]
-         -D INT     max total depth [$opts{D}]
-         -q INT     min SNP quality [$opts{q}]
-         -Q INT     min RMS mapQ (given MQ) [$opts{Q}]
-         -1 FLOAT   min P-value for strand bias (given PV4) [$opts{1}]
-         -2 FLOAT   min P-value for baseQ bias [$opts{2}]
-         -3 FLOAT   min P-value for mapQ bias [$opts{3}]
-         -4 FLOAT   min P-value for end distance bias [$opts{4}]\n
-/) if (@ARGV == 0 && -t STDIN);
-
-  my %ts = (AG=>1, GA=>1, CT=>1, TC=>1);
-
-  my @n = (0, 0);
-  while (<>) {
-       if (/^#/) {
-         print;
-         next;
-       }
-       next if (/PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/ && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
-       my $depth = -1;
-       $depth = $1 if (/DP=(\d+)/);
-       $depth = $1+$2+$3+$4 if (/DP4=(\d+),(\d+),(\d+),(\d+)/);
-       next if ($depth > 0 && ($depth < $opts{d} || $depth > $opts{D}));
-       next if (/MQ=(\d+)/ && $1 < $opts{Q});
-       my @t = split;
-       next if ($t[5] >= 0 && $t[5] < $opts{q});
-       ++$n[0];
-       my @s = split(',', $t[4]);
-       ++$n[1] if ($ts{$t[3].$s[0]});
-       print;
+       print STDERR join("\t", substr("UQdDaGgP", $first->[1], 1), @$first[3 .. @$first-1]), "\n";

   }
 }
 
   }
 }
 


@@ -436,7 +406,6 @@ Command: subsam       get a subset of samples
          fillac       fill the allele count field
          qstats       SNP stats stratified by QUAL
          varFilter    filtering short variants
          fillac       fill the allele count field
          qstats       SNP stats stratified by QUAL
          varFilter    filtering short variants

-         filter4vcf   filtering VCFs produced by samtools+bcftools

          hapmap2vcf   convert the hapmap format to VCF
          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 \n/);
          hapmap2vcf   convert the hapmap format to VCF
          ucscsnp2vcf  convert UCSC SNP SQL dump to VCF
 \n/);