X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fcall1.c;h=20083016e2814e51dd0fda5e8e4caa78248fe8bd;hb=28097a71d7a946ceb6e2091b724a7572b5786c95;hp=22ff2aca91be76f564e7bd3825abddcd4616ba42;hpb=faf76ae18e5ddf2d03253e2ed77f444bcdf4914c;p=samtools.git

diff --git a/bcftools/call1.c b/bcftools/call1.c
index 22ff2ac..2008301 100644
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -190,7 +190,27 @@ static char **read_samples(const char *fn, int *_n)
 	*_n = 0;
 	s.l = s.m = 0; s.s = 0;
 	fp = gzopen(fn, "r");
-	if (fp == 0) return 0; // fail to open file
+	if (fp == 0) 
+    {
+        // interpret as sample names, not as a file name
+        const char *t = fn, *p = t;
+        while (*t)
+        {
+            t++;
+            if ( *t==',' || !*t )
+            { 
+                sam = realloc(sam, sizeof(void*)*(n+1));
+                sam[n] = (char*) malloc(sizeof(char)*(t-p+2));
+                memcpy(sam[n], p, t-p);
+                sam[n][t-p]   = 0;
+                sam[n][t-p+1] = 2;    // assume diploid
+                p = t+1;
+                n++; 
+            }
+        }
+        *_n = n;
+        return sam; // fail to open file
+    }
 	ks = ks_init(fp);
 	while (ks_getuntil(ks, 0, &s, &dret) >= 0) {
 		int l;
@@ -238,8 +258,10 @@ static void write_header(bcf_hdr_t *h)
 		kputs("##INFO=<ID=AC1,Number=1,Type=Float,Description=\"Max-likelihood estimate of the first ALT allele count (no HWE assumption)\">\n", &str);
 	if (!strstr(str.s, "##INFO=<ID=AN,"))
 		kputs("##INFO=<ID=AN,Number=1,Type=Integer,Description=\"Total number of alleles in called genotypes\">\n", &str);
-	if (!strstr(str.s, "##INFO=<ID=IS,"))
-		kputs("##INFO=<ID=IS,Number=2,Type=Float,Description=\"Maximum number of reads supporting an indel and fraction of indel reads\">\n", &str);
+	if (!strstr(str.s, "##INFO=<ID=IDV,"))
+		kputs("##INFO=<ID=IDV,Number=1,Type=Integer,Description=\"Maximum number of reads supporting an indel\">\n", &str);
+	if (!strstr(str.s, "##INFO=<ID=IMF,"))
+		kputs("##INFO=<ID=IMF,Number=1,Type=Float,Description=\"Maximum fraction of reads supporting an indel\">\n", &str);
 	if (!strstr(str.s, "##INFO=<ID=AC,"))
 		kputs("##INFO=<ID=AC,Number=A,Type=Integer,Description=\"Allele count in genotypes for each ALT allele, in the same order as listed\">\n", &str);
 	if (!strstr(str.s, "##INFO=<ID=G3,"))
@@ -266,8 +288,16 @@ 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=QBD,"))
+        kputs("##INFO=<ID=QBD,Number=1,Type=Float,Description=\"Quality by Depth: QUAL/#reads\">\n", &str);
+    //if (!strstr(str.s, "##INFO=<ID=RPS,"))
+    //    kputs("##INFO=<ID=RPS,Number=3,Type=Float,Description=\"Read Position Stats: depth, average, stddev\">\n", &str);
+    if (!strstr(str.s, "##INFO=<ID=RPB,"))
+        kputs("##INFO=<ID=RPB,Number=1,Type=Float,Description=\"Read Position Bias\">\n", &str);
+    if (!strstr(str.s, "##INFO=<ID=MDV,"))
+        kputs("##INFO=<ID=MDV,Number=1,Type=Integer,Description=\"Maximum number of high-quality nonRef reads in samples\">\n", &str);
     if (!strstr(str.s, "##INFO=<ID=VDB,"))
-        kputs("##INFO=<ID=VDB,Number=1,Type=Float,Description=\"Variant Distance Bias\">\n", &str);
+        kputs("##INFO=<ID=VDB,Number=1,Type=Float,Description=\"Variant Distance Bias (v2) for filtering splice-site artefacts in RNA-seq data. Note: this version may be broken.\">\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,"))
@@ -298,6 +328,9 @@ int bcfview(int argc, char *argv[])
 	extern int bcf_trio_call(uint32_t *prep, const bcf1_t *b, int *llr, int64_t *gt);
 	extern int bcf_pair_call(const bcf1_t *b);
 	extern int bcf_min_diff(const bcf1_t *b);
+	extern int bcf_p1_get_M(bcf_p1aux_t *b);
+
+	extern gzFile bcf_p1_fp_lk;
 
 	bcf_t *bp, *bout = 0;
 	bcf1_t *b, *blast;
@@ -313,10 +346,10 @@ 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; vc.min_ma_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:Ywm:")) >= 0) {
+	while ((c = getopt(argc, argv, "FN1:l:cC:eHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:T:Ywm:K:")) >= 0) {
 		switch (c) {
 		case '1': vc.n1 = atoi(optarg); break;
-		case 'l': vc.bed = bed_read(optarg); break;
+		case 'l': vc.bed = bed_read(optarg); if (!vc.bed) { fprintf(stderr,"Could not read \"%s\"\n", optarg); return 1; } break;
 		case 'D': vc.fn_dict = strdup(optarg); break;
 		case 'F': vc.flag |= VC_FIX_PL; break;
 		case 'N': vc.flag |= VC_ACGT_ONLY; break;
@@ -343,6 +376,7 @@ int bcfview(int argc, char *argv[])
 		case 'C': vc.min_lrt = atof(optarg); break;
 		case 'X': vc.min_perm_p = atof(optarg); break;
 		case 'd': vc.min_smpl_frac = atof(optarg); break;
+		case 'K': bcf_p1_fp_lk = gzopen(optarg, "w"); 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];
@@ -432,7 +466,7 @@ int bcfview(int argc, char *argv[])
 			vc.sublist = calloc(vc.n_sub, sizeof(int));
 			hout = bcf_hdr_subsam(hin, vc.n_sub, vc.subsam, vc.sublist);
 		}
-		if (vc.flag & VC_CALL) write_header(hout);
+		write_header(hout); // always print the header
 		vcf_hdr_write(bout, hout);
 	}
 	if (vc.flag & VC_CALL) {
@@ -466,6 +500,10 @@ int bcfview(int argc, char *argv[])
 			}
 		}
 	}
+	if (bcf_p1_fp_lk && p1) {
+		int32_t M = bcf_p1_get_M(p1);
+		gzwrite(bcf_p1_fp_lk, &M, 4);
+	}
 	while (vcf_read(bp, hin, b) > 0) {
 		int is_indel, cons_llr = -1;
 		int64_t cons_gt = -1;
@@ -514,15 +552,19 @@ int bcfview(int argc, char *argv[])
 			int i;
 			for (i = 0; i < 9; ++i) em[i] = -1.;
 		}
-        if ( !(vc.flag&VC_KEEPALT) && vc.flag&VC_CALL && vc.min_ma_lrt>=0 )
+        if ( !(vc.flag&VC_KEEPALT) && (vc.flag&VC_CALL) && vc.min_ma_lrt>=0 )
         {
             bcf_p1_set_ploidy(b, p1); // could be improved: do this per site to allow pseudo-autosomal regions
-            int gts = call_multiallelic_gt(b,p1,vc.min_ma_lrt);
+            int gts = call_multiallelic_gt(b, p1, vc.min_ma_lrt, vc.flag&VC_VARONLY);
             if ( gts<=1 && vc.flag & VC_VARONLY ) continue;
         }
 		else if (vc.flag & VC_CALL) { // call variants
 			bcf_p1rst_t pr;
-			int calret = bcf_p1_cal(b, (em[7] >= 0 && em[7] < vc.min_lrt), p1, &pr);
+			int calret;
+			gzwrite(bcf_p1_fp_lk, &b->tid, 4);
+			gzwrite(bcf_p1_fp_lk, &b->pos, 4);
+			gzwrite(bcf_p1_fp_lk, &em[0], sizeof(double));
+			calret = bcf_p1_cal(b, (em[7] >= 0 && em[7] < vc.min_lrt), p1, &pr);
 			if (n_processed % 100000 == 0) {
 				fprintf(stderr, "[%s] %ld sites processed.\n", __func__, (long)n_processed);
 				bcf_p1_dump_afs(p1);
@@ -567,6 +609,8 @@ int bcfview(int argc, char *argv[])
 		} else bcf_fix_gt(b);
 		vcf_write(bout, hout, b);
 	}
+
+	if (bcf_p1_fp_lk) gzclose(bcf_p1_fp_lk);
 	if (vc.prior_file) free(vc.prior_file);
 	if (vc.flag & VC_CALL) bcf_p1_dump_afs(p1);
 	if (hin != hout) bcf_hdr_destroy(hout);