X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fcall1.c;h=6c530085ba591a23aeae98918361cb10745793a6;hb=ccb7838cb53bf0ca6917a77f7991d940057c12db;hp=b0baee2ea1b346e24eadd2d9bccea0ed6e886f6a;hpb=23e43ef6778e51eb2a71046282d1c21431520857;p=samtools.git

diff --git a/bcftools/call1.c b/bcftools/call1.c
index b0baee2..6c53008 100644
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -33,6 +33,7 @@ KSTREAM_INIT(gzFile, gzread, 16384)
 #define VC_EM       0x10000
 #define VC_PAIRCALL 0x20000
 #define VC_QCNT     0x40000
+#define VC_INDEL_ONLY 0x80000
 
 typedef struct {
 	int flag, prior_type, n1, n_sub, *sublist, n_perm;
@@ -131,7 +132,6 @@ static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr,
 		if (em[5] >= 0 && em[6] >= 0) ksprintf(&s, ";AF2=%.4g,%.4g", 1 - em[5], 1 - em[6]);
 		if (em[7] >= 0) ksprintf(&s, ";LRT=%.3g", em[7]);
 		if (em[8] >= 0) ksprintf(&s, ";LRT2=%.3g", em[8]);
-		//if (em[9] >= 0) ksprintf(&s, ";LRT1=%.3g", em[9]);
 	}
 	if (cons_llr > 0) {
 		ksprintf(&s, ";CLR=%d", cons_llr);
@@ -274,6 +274,8 @@ static void write_header(bcf_hdr_t *h)
         kputs("##INFO=<ID=QCHI2,Number=1,Type=Integer,Description=\"Phred scaled PCHI2.\">\n", &str);
     if (!strstr(str.s, "##INFO=<ID=RP,"))
         kputs("##INFO=<ID=PR,Number=1,Type=Integer,Description=\"# permutations yielding a smaller PCHI2.\">\n", &str);
+    if (!strstr(str.s, "##INFO=<ID=VDB,"))
+        kputs("##INFO=<ID=VDB,Number=1,Type=Float,Description=\"Variant Distance Bias\">\n", &str);
     if (!strstr(str.s, "##FORMAT=<ID=GT,"))
         kputs("##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n", &str);
     if (!strstr(str.s, "##FORMAT=<ID=GQ,"))
@@ -282,10 +284,12 @@ static void write_header(bcf_hdr_t *h)
         kputs("##FORMAT=<ID=GL,Number=3,Type=Float,Description=\"Likelihoods for RR,RA,AA genotypes (R=ref,A=alt)\">\n", &str);
 	if (!strstr(str.s, "##FORMAT=<ID=DP,"))
 		kputs("##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"# high-quality bases\">\n", &str);
+	if (!strstr(str.s, "##FORMAT=<ID=DV,"))
+		kputs("##FORMAT=<ID=DV,Number=1,Type=Integer,Description=\"# high-quality non-reference bases\">\n", &str);
 	if (!strstr(str.s, "##FORMAT=<ID=SP,"))
 		kputs("##FORMAT=<ID=SP,Number=1,Type=Integer,Description=\"Phred-scaled strand bias P-value\">\n", &str);
 	if (!strstr(str.s, "##FORMAT=<ID=PL,"))
-		kputs("##FORMAT=<ID=PL,Number=-1,Type=Integer,Description=\"List of Phred-scaled genotype likelihoods, number of values is (#ALT+1)*(#ALT+2)/2\">\n", &str);
+		kputs("##FORMAT=<ID=PL,Number=G,Type=Integer,Description=\"List of Phred-scaled genotype likelihoods\">\n", &str);
 	h->l_txt = str.l + 1; h->txt = str.s;
 }
 
@@ -317,7 +321,7 @@ int bcfview(int argc, char *argv[])
 	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; vc.min_smpl_frac = 0; vc.min_lrt = 1;
 	memset(qcnt, 0, 8 * 256);
-	while ((c = getopt(argc, argv, "FN1:l:cC:eHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:T:Y")) >= 0) {
+	while ((c = getopt(argc, argv, "FN1:l:cC:eHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:T:Yw")) >= 0) {
 		switch (c) {
 		case '1': vc.n1 = atoi(optarg); break;
 		case 'l': vc.bed = bed_read(optarg); break;
@@ -334,6 +338,7 @@ int bcfview(int argc, char *argv[])
 		case 'u': vc.flag |= VC_UNCOMP | VC_BCFOUT; break;
 		case 'g': vc.flag |= VC_CALL_GT | VC_CALL; break;
 		case 'I': vc.flag |= VC_NO_INDEL; break;
+		case 'w': vc.flag |= VC_INDEL_ONLY; break;
 		case 'M': vc.flag |= VC_ANNO_MAX; break;
 		case 'Y': vc.flag |= VC_QCNT; break;
 		case 't': vc.theta = atof(optarg); break;
@@ -481,6 +486,7 @@ int bcfview(int argc, char *argv[])
 		if (vc.flag & VC_FIX_PL) bcf_fix_pl(b);
 		is_indel = bcf_is_indel(b);
 		if ((vc.flag & VC_NO_INDEL) && is_indel) continue;
+		if ((vc.flag & VC_INDEL_ONLY) && !is_indel) continue;
 		if ((vc.flag & VC_ACGT_ONLY) && !is_indel) {
 			int x;
 			if (b->ref[0] == 0 || b->ref[1] != 0) continue;