X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fcall1.c;h=b2d7d4a197b112ea10c6caf3ad4a9a13d95b1298;hb=3384f6c6c1642ada4edae9204ca1202672de7d5a;hp=1c35ee85a64bcf607d4f7b8678280329cb712a9d;hpb=4d48006473298b62ee54f13d706f9f2e512e02bd;p=samtools.git

diff --git a/bcftools/call1.c b/bcftools/call1.c
index 1c35ee8..b2d7d4a 100644
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -26,9 +26,12 @@ KSTREAM_INIT(gzFile, gzread, 16384)
 #define VC_ANNO_MAX 16384
 #define VC_FIX_PL   32768
 #define VC_EM       0x10000
+#define VC_PAIRCALL 0x20000
+#define VC_QCNT     0x40000
 
 typedef struct {
 	int flag, prior_type, n1, n_sub, *sublist, n_perm;
+	uint32_t *trio_aux;
 	char *prior_file, **subsam, *fn_dict;
 	uint8_t *ploidy;
 	double theta, pref, indel_frac, min_perm_p, min_smpl_frac, min_lrt;
@@ -102,7 +105,7 @@ static void rm_info(bcf1_t *b, const char *key)
 	bcf_sync(b);
 }
 
-static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, double pref, int flag, double em[10])
+static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, double pref, int flag, double em[10], int cons_llr, int64_t cons_gt)
 {
 	kstring_t s;
 	int has_I16, is_var;
@@ -125,6 +128,12 @@ static int update_bcf1(bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr,
 		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);
+		if (cons_gt > 0)
+			ksprintf(&s, ";UGT=%c%c%c;CGT=%c%c%c", cons_gt&0xff, cons_gt>>8&0xff, cons_gt>>16&0xff,
+				     cons_gt>>32&0xff, cons_gt>>40&0xff, cons_gt>>48&0xff);
+	}
 	if (pr == 0) { // if pr is unset, return
 		kputc('\0', &s); kputs(b->fmt, &s); kputc('\0', &s);
 		free(b->str);
@@ -240,6 +249,12 @@ static void write_header(bcf_hdr_t *h)
 		kputs("##INFO=<ID=G3,Number=3,Type=Float,Description=\"ML estimate of genotype frequencies\">\n", &str);
 	if (!strstr(str.s, "##INFO=<ID=HWE,"))
 		kputs("##INFO=<ID=HWE,Number=1,Type=Float,Description=\"Chi^2 based HWE test P-value based on G3\">\n", &str);
+	if (!strstr(str.s, "##INFO=<ID=CLR,"))
+		kputs("##INFO=<ID=CLR,Number=1,Type=Integer,Description=\"Log ratio of genotype likelihoods with and without the constraint\">\n", &str);
+	if (!strstr(str.s, "##INFO=<ID=UGT,"))
+		kputs("##INFO=<ID=UGT,Number=1,Type=String,Description=\"The most probable unconstrained genotype configuration in the trio\">\n", &str);
+	if (!strstr(str.s, "##INFO=<ID=CGT,"))
+		kputs("##INFO=<ID=CGT,Number=1,Type=String,Description=\"The most probable constrained genotype configuration in the trio\">\n", &str);
 //	if (!strstr(str.s, "##INFO=<ID=CI95,"))
 //		kputs("##INFO=<ID=CI95,Number=2,Type=Float,Description=\"Equal-tail Bayesian credible interval of the site allele frequency at the 95% level\">\n", &str);
 	if (!strstr(str.s, "##INFO=<ID=PV4,"))
@@ -278,10 +293,15 @@ int bcfview(int argc, char *argv[])
 	extern int bcf_fix_gt(bcf1_t *b);
 	extern int bcf_anno_max(bcf1_t *b);
 	extern int bcf_shuffle(bcf1_t *b, int seed);
+	extern uint32_t *bcf_trio_prep(int is_x, int is_son);
+	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);
+
 	bcf_t *bp, *bout = 0;
 	bcf1_t *b, *blast;
 	int c, *seeds = 0;
-	uint64_t n_processed = 0;
+	uint64_t n_processed = 0, qcnt[256];
 	viewconf_t vc;
 	bcf_p1aux_t *p1 = 0;
 	bcf_hdr_t *hin, *hout;
@@ -291,7 +311,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; vc.min_smpl_frac = 0; vc.min_lrt = 1;
-	while ((c = getopt(argc, argv, "FN1:l:cC:eHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:")) >= 0) {
+	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) {
 		switch (c) {
 		case '1': vc.n1 = atoi(optarg); break;
 		case 'l': vc.bed = bed_read(optarg); break;
@@ -309,6 +330,7 @@ int bcfview(int argc, char *argv[])
 		case 'g': vc.flag |= VC_CALL_GT | VC_CALL; break;
 		case 'I': vc.flag |= VC_NO_INDEL; break;
 		case 'M': vc.flag |= VC_ANNO_MAX; break;
+		case 'Y': vc.flag |= VC_QCNT; break;
 		case 't': vc.theta = atof(optarg); break;
 		case 'p': vc.pref = atof(optarg); break;
 		case 'i': vc.indel_frac = atof(optarg); break;
@@ -323,6 +345,16 @@ int bcfview(int argc, char *argv[])
 			for (tid = 0; tid < vc.n_sub; ++tid) vc.ploidy[tid] = vc.subsam[tid][strlen(vc.subsam[tid]) + 1];
 			tid = -1;
 			break;
+		case 'T':
+			if (strcmp(optarg, "trioauto") == 0) vc.trio_aux = bcf_trio_prep(0, 0);
+			else if (strcmp(optarg, "trioxd") == 0) vc.trio_aux = bcf_trio_prep(1, 0);
+			else if (strcmp(optarg, "trioxs") == 0) vc.trio_aux = bcf_trio_prep(1, 1);
+			else if (strcmp(optarg, "pair") == 0) vc.flag |= VC_PAIRCALL;
+			else {
+				fprintf(stderr, "[%s] Option '-T' can only take value trioauto, trioxd or trioxs.\n", __func__);
+				return 1;
+			}
+			break;
 		case 'P':
 			if (strcmp(optarg, "full") == 0) vc.prior_type = MC_PTYPE_FULL;
 			else if (strcmp(optarg, "cond2") == 0) vc.prior_type = MC_PTYPE_COND2;
@@ -357,6 +389,7 @@ int bcfview(int argc, char *argv[])
 		fprintf(stderr, "       -p FLOAT  variant if P(ref|D)<FLOAT [%.3g]\n", vc.pref);
 		fprintf(stderr, "       -P STR    type of prior: full, cond2, flat [full]\n");
 		fprintf(stderr, "       -t FLOAT  scaled substitution mutation rate [%.4g]\n", vc.theta);
+		fprintf(stderr, "       -T STR    constrained calling; STR can be: pair, trioauto, trioxd and trioxs (see manual) [null]\n");
 		fprintf(stderr, "       -v        output potential variant sites only (force -c)\n");
 		fprintf(stderr, "\nContrast calling and association test options:\n\n");
 		fprintf(stderr, "       -1 INT    number of group-1 samples [0]\n");
@@ -430,7 +463,8 @@ int bcfview(int argc, char *argv[])
 		}
 	}
 	while (vcf_read(bp, hin, b) > 0) {
-		int is_indel;
+		int is_indel, cons_llr = -1;
+		int64_t cons_gt = -1;
 		double em[10];
 		if ((vc.flag & VC_VARONLY) && strcmp(b->alt, "X") == 0) continue;
 		if ((vc.flag & VC_VARONLY) && vc.min_smpl_frac > 0.) {
@@ -456,10 +490,19 @@ int bcfview(int argc, char *argv[])
 			if (!(l > begin && end > b->pos)) continue;
 		}
 		++n_processed;
+		if (vc.flag & VC_QCNT) { // summarize the difference
+			int x = bcf_min_diff(b);
+			if (x > 255) x = 255;
+			if (x >= 0) ++qcnt[x];
+		}
 		if (vc.flag & VC_QCALL) { // output QCALL format; STOP here
 			bcf_2qcall(hout, b);
 			continue;
 		}
+		if (vc.trio_aux) // do trio calling
+			bcf_trio_call(vc.trio_aux, b, &cons_llr, &cons_gt);
+		else if (vc.flag & VC_PAIRCALL)
+			cons_llr = bcf_pair_call(b);
 		if (vc.flag & (VC_CALL|VC_ADJLD|VC_EM)) bcf_gl2pl(b);
 		if (vc.flag & VC_EM) bcf_em1(b, vc.n1, 0x1ff, em);
 		else {
@@ -491,8 +534,8 @@ int bcfview(int argc, char *argv[])
 				}
 				pr.perm_rank = n;
 			}
-			if (calret >= 0) update_bcf1(b, p1, &pr, vc.pref, vc.flag, em);
-		} else if (vc.flag & VC_EM) update_bcf1(b, 0, 0, 0, vc.flag, em);
+			if (calret >= 0) update_bcf1(b, p1, &pr, vc.pref, vc.flag, em, cons_llr, cons_gt);
+		} else if (vc.flag & VC_EM) update_bcf1(b, 0, 0, 0, vc.flag, em, cons_llr, cons_gt);
 		if (vc.flag & VC_ADJLD) { // compute LD
 			double f[4], r2;
 			if ((r2 = bcf_pair_freq(blast, b, f)) >= 0) {
@@ -521,12 +564,16 @@ int bcfview(int argc, char *argv[])
 	vcf_close(bp); vcf_close(bout);
 	if (vc.fn_dict) free(vc.fn_dict);
 	if (vc.ploidy) free(vc.ploidy);
+	if (vc.trio_aux) free(vc.trio_aux);
 	if (vc.n_sub) {
 		int i;
 		for (i = 0; i < vc.n_sub; ++i) free(vc.subsam[i]);
 		free(vc.subsam); free(vc.sublist);
 	}
 	if (vc.bed) bed_destroy(vc.bed);
+	if (vc.flag & VC_QCNT)
+		for (c = 0; c < 256; ++c)
+			fprintf(stderr, "QT\t%d\t%lld\n", c, (long long)qcnt[c]);
 	if (seeds) free(seeds);
 	if (p1) bcf_p1_destroy(p1);
 	return 0;