Release samtools-0.1.14 (r933:170)

[samtools.git] / bcftools / call1.c

diff --git a/bcftools/call1.c b/bcftools/call1.c
index a520e3cc966553c2acfbef98a0ddfc2e4c48b53f..e074fb5cac2fec5bf2129aec0c2c736764dd0ef5 100644 (file)
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -33,7 +33,7 @@ typedef struct {
        int flag, prior_type, n1, n_sub, *sublist, n_perm;
        char *fn_list, *prior_file, **subsam, *fn_dict;
        uint8_t *ploidy;
-       double theta, pref, indel_frac, min_perm_p;
+       double theta, pref, indel_frac, min_perm_p, min_smpl_frac;
 } viewconf_t;
 
 khash_t(set64) *bcf_load_pos(const char *fn, bcf_hdr_t *_h)
@@ -297,8 +297,8 @@ int bcfview(int argc, char *argv[])
 
        tid = begin = end = -1;
        memset(&vc, 0, sizeof(viewconf_t));
-       vc.prior_type = vc.n1 = -1; vc.theta = 1e-3; vc.pref = 0.5; vc.indel_frac = -1.; vc.n_perm = 0; vc.min_perm_p = 0.01;
-       while ((c = getopt(argc, argv, "FN1:l:cHAGvbSuP:t:p:QgLi:IMs:D:U:X:")) >= 0) {
+       vc.prior_type = vc.n1 = -1; vc.theta = 1e-3; vc.pref = 0.5; vc.indel_frac = -1.; vc.n_perm = 0; vc.min_perm_p = 0.01; vc.min_smpl_frac = 0;
+       while ((c = getopt(argc, argv, "FN1:l:cHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:")) >= 0) {
                switch (c) {
                case '1': vc.n1 = atoi(optarg); break;
                case 'l': vc.fn_list = strdup(optarg); break;
@@ -322,6 +322,7 @@ int bcfview(int argc, char *argv[])
                case 'L': vc.flag |= VC_ADJLD; break;
                case 'U': vc.n_perm = atoi(optarg); break;
                case 'X': vc.min_perm_p = atof(optarg); break;
+               case 'd': vc.min_smpl_frac = atof(optarg); break;
                case 's': vc.subsam = read_samples(optarg, &vc.n_sub);
                        vc.ploidy = calloc(vc.n_sub + 1, 1);
                        for (tid = 0; tid < vc.n_sub; ++tid) vc.ploidy[tid] = vc.subsam[tid][strlen(vc.subsam[tid]) + 1];
@@ -351,6 +352,7 @@ int bcfview(int argc, char *argv[])
                fprintf(stderr, "         -L        calculate LD for adjacent sites\n");
                fprintf(stderr, "         -I        skip indels\n");
                fprintf(stderr, "         -F        PL generated by r921 or before\n");
+               fprintf(stderr, "         -d FLOAT  skip loci where less than FLOAT fraction of samples covered [0]\n");
                fprintf(stderr, "         -D FILE   sequence dictionary for VCF->BCF conversion [null]\n");
                fprintf(stderr, "         -l FILE   list of sites to output [all sites]\n");
                fprintf(stderr, "         -t FLOAT  scaled substitution mutation rate [%.4g]\n", vc.theta);
@@ -429,6 +431,12 @@ int bcfview(int argc, char *argv[])
        }
        while (vcf_read(bp, hin, b) > 0) {
                int is_indel;
+               if ((vc.flag & VC_VARONLY) && strcmp(b->alt, "X") == 0) continue;
+               if ((vc.flag & VC_VARONLY) && vc.min_smpl_frac > 0.) {
+                       extern int bcf_smpl_covered(const bcf1_t *b);
+                       int n = bcf_smpl_covered(b);
+                       if ((double)n / b->n_smpl < vc.min_smpl_frac) continue;
+               }
                if (vc.n_sub) bcf_subsam(vc.n_sub, vc.sublist, b);
                if (vc.flag & VC_FIX_PL) bcf_fix_pl(b);
                is_indel = bcf_is_indel(b);