return 0;
}
-int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b)
+int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bcr)
{
+ extern double kt_fisher_exact(int n11, int n12, int n21, int n22, double *_left, double *_right, double *two);
kstring_t s;
int i;
b->n_smpl = bc->n;
}
kputc('\0', &s);
// FMT
- kputs("PL", &s); kputc('\0', &s);
+ kputs("PL", &s);
+ if (bcr) {
+ kputs(":DP", &s);
+ }
+ kputc('\0', &s);
b->m_str = s.m; b->str = s.s; b->l_str = s.l;
bcf_sync(b);
memcpy(b->gi[0].data, bc->PL, b->gi[0].len * bc->n);
+ if (bcr) {
+ uint16_t *dp = (uint16_t*)b->gi[1].data;
+ for (i = 0; i < bc->n; ++i) {
+ dp[i] = bcr[i].depth < 0xffff? bcr[i].depth : 0xffff;
+ }
+ }
return 0;
}
void bcf_call_destroy(bcf_callaux_t *bca);
int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base /*4-bit*/, bcf_callaux_t *bca, bcf_callret1_t *r);
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);
+ int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bcr);
#ifdef __cplusplus
}
#define MPLP_NO_COMP 0x20
#define MPLP_NO_ORPHAN 0x40
#define MPLP_REALN 0x80
+#define MPLP_SAMINFO 0x100
typedef struct {
int max_mq, min_mq, flag, min_baseQ, capQ_thres;
for (i = 0; i < gplp.n; ++i)
bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], ref16, bca, bcr + i);
bcf_call_combine(gplp.n, bcr, ref16, &bc);
- bcf_call2bcf(tid, pos, &bc, b);
+ bcf_call2bcf(tid, pos, &bc, b, (conf->flag&MPLP_SAMINFO)? bcr : 0);
bcf_write(bp, bh, b);
bcf_destroy(b);
} else {
mplp.min_baseQ = 13;
mplp.capQ_thres = 0;
mplp.flag = MPLP_NO_ORPHAN | MPLP_REALN;
- while ((c = getopt(argc, argv, "gf:r:l:M:q:Q:uaORC:B")) >= 0) {
+ while ((c = getopt(argc, argv, "gf:r:l:M:q:Q:uaORC:BI")) >= 0) {
switch (c) {
case 'f':
mplp.fai = fai_load(optarg);
case 'B': mplp.flag &= ~MPLP_REALN & ~MPLP_NO_ORPHAN; break;
case 'O': mplp.flag |= MPLP_NO_ORPHAN; break;
case 'R': mplp.flag |= MPLP_REALN; break;
+ case 'I': mplp.flag |= MPLP_SAMINFO; break;
case 'C': mplp.capQ_thres = atoi(optarg); break;
case 'M': mplp.max_mq = atoi(optarg); break;
case 'q': mplp.min_mq = atoi(optarg); break;
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, " -I additional individual information (DP and ST)\n");
fprintf(stderr, "\n");
fprintf(stderr, "Notes: Assuming diploid individuals.\n\n");
return 1;
b->gi[i].len = b->n_alleles * (b->n_alleles + 1) / 2;
} else if (b->gi[i].fmt == bcf_str2int("DP", 2) || b->gi[i].fmt == bcf_str2int("HQ", 2)) {
b->gi[i].len = 2;
- } else if (b->gi[i].fmt == bcf_str2int("GQ", 2) || b->gi[i].fmt == bcf_str2int("GT", 2)) {
+ } else if (b->gi[i].fmt == bcf_str2int("GQ", 2) || b->gi[i].fmt == bcf_str2int("GT", 2)
+ || b->gi[i].fmt == bcf_str2int("ST", 2))
+ {
b->gi[i].len = 1;
} else if (b->gi[i].fmt == bcf_str2int("GL", 2)) {
b->gi[i].len = b->n_alleles * (b->n_alleles + 1) / 2 * 4;
}
} else if (b->gi[i].fmt == bcf_str2int("DP", 2)) {
kputw(((uint16_t*)b->gi[i].data)[j], s);
- } else if (b->gi[i].fmt == bcf_str2int("GQ", 2)) {
+ } else if (b->gi[i].fmt == bcf_str2int("GQ", 2) || b->gi[i].fmt == bcf_str2int("ST", 2)) {
kputw(((uint8_t*)b->gi[i].data)[j], s);
} else if (b->gi[i].fmt == bcf_str2int("GT", 2)) {
int y = ((uint8_t*)b->gi[i].data)[j];