From: Heng Li Date: Tue, 12 Oct 2010 19:48:46 +0000 (+0000) Subject: * vcfutils.pl: fixed a typo in help message X-Git-Url: https://git.donarmstrong.com/?p=samtools.git;a=commitdiff_plain;h=032937695157a15fb0ce042deda22afc30be1616 * vcfutils.pl: fixed a typo in help message * added APIs: bcf_append_info() and bcf_cpy() * calculate adjacent LD --- diff --git a/bcftools/Makefile b/bcftools/Makefile index 2dbcf59..b6ab972 100644 --- a/bcftools/Makefile +++ b/bcftools/Makefile @@ -1,7 +1,7 @@ CC= gcc -CFLAGS= -g -Wall -O2 #-m64 #-arch ppc +CFLAGS= -g -Wall #-O2 #-m64 #-arch ppc DFLAGS= -D_FILE_OFFSET_BITS=64 -D_USE_KNETFILE -LOBJS= bcf.o vcf.o bcfutils.o prob1.o kfunc.o index.o fet.o bcf2qcall.o +LOBJS= bcf.o vcf.o bcfutils.o prob1.o ld.o kfunc.o index.o fet.o bcf2qcall.o OMISC= .. AOBJS= call1.o main.o $(OMISC)/kstring.o $(OMISC)/bgzf.o $(OMISC)/knetfile.o PROG= bcftools diff --git a/bcftools/bcf.c b/bcftools/bcf.c index 35c3d94..3531647 100644 --- a/bcftools/bcf.c +++ b/bcftools/bcf.c @@ -267,3 +267,41 @@ char *bcf_fmt(const bcf_hdr_t *h, bcf1_t *b) bcf_fmt_core(h, b, &s); return s.s; } + +int bcf_append_info(bcf1_t *b, const char *info, int l) +{ + int shift = b->fmt - b->str; + int l_fmt = b->l_str - shift; + char *ori = b->str; + if (b->l_str + l > b->m_str) { // enlarge if necessary + b->m_str = b->l_str + l; + kroundup32(b->m_str); + b->str = realloc(b->str, b->m_str); + } + memmove(b->str + shift + l, b->str + shift, l_fmt); // move the FORMAT field + memcpy(b->str + shift - 1, info, l); // append to the INFO field + b->str[shift + l - 1] = '\0'; + b->fmt = b->str + shift + l; + b->l_str += l; + bcf_sync(b); +// if (ori != b->str) bcf_sync(b); // synchronize when realloc changes the pointer + return 0; +} + +int bcf_cpy(bcf1_t *r, const bcf1_t *b) +{ + char *t1 = r->str; + bcf_ginfo_t *t2 = r->gi; + int i, t3 = r->m_str, t4 = r->m_gi; + *r = *b; + r->str = t1; r->gi = t2; r->m_str = t3; r->m_gi = t4; + if (r->m_str < b->m_str) { + r->m_str = b->m_str; + r->str = realloc(r->str, r->m_str); + } + memcpy(r->str, b->str, r->m_str); + bcf_sync(r); // calling bcf_sync() is simple but inefficient + for (i = 0; i < r->n_gi; ++i) + memcpy(r->gi[i].data, b->gi[i].data, r->n_smpl * r->gi[i].len); + return 0; +} diff --git a/bcftools/bcf.h b/bcftools/bcf.h index 7286111..3657895 100644 --- a/bcftools/bcf.h +++ b/bcftools/bcf.h @@ -106,6 +106,10 @@ extern "C" { int bcf_destroy(bcf1_t *b); // BCF->VCF conversion char *bcf_fmt(const bcf_hdr_t *h, bcf1_t *b); + // append more info + int bcf_append_info(bcf1_t *b, const char *info, int l); + // copy + int bcf_cpy(bcf1_t *r, const bcf1_t *b); // open a VCF or BCF file if "b" is set in "mode" bcf_t *vcf_open(const char *fn, const char *mode); diff --git a/bcftools/call1.c b/bcftools/call1.c index 1eacd4f..aa3e591 100644 --- a/bcftools/call1.c +++ b/bcftools/call1.c @@ -24,6 +24,7 @@ KSTREAM_INIT(gzFile, gzread, 16384) #define VC_ACGT_ONLY 512 #define VC_QCALL 1024 #define VC_CALL_GT 2048 +#define VC_ADJLD 4096 typedef struct { int flag, prior_type, n1; @@ -193,11 +194,13 @@ static int update_bcf1(int n_smpl, bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p return is_var; } +double bcf_ld_freq(const bcf1_t *b0, const bcf1_t *b1, double f[4]); + int bcfview(int argc, char *argv[]) { extern int bcf_2qcall(bcf_hdr_t *h, bcf1_t *b); bcf_t *bp, *bout = 0; - bcf1_t *b; + bcf1_t *b, *blast; int c; uint64_t n_processed = 0; viewconf_t vc; @@ -210,7 +213,7 @@ 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; - while ((c = getopt(argc, argv, "N1:l:cHAGvbSuP:t:p:Qg")) >= 0) { + while ((c = getopt(argc, argv, "N1:l:cHAGvbSuP:t:p:QgL")) >= 0) { switch (c) { case '1': vc.n1 = atoi(optarg); break; case 'l': vc.fn_list = strdup(optarg); break; @@ -227,6 +230,7 @@ int bcfview(int argc, char *argv[]) case 't': vc.theta = atof(optarg); break; case 'p': vc.pref = atof(optarg); break; case 'Q': vc.flag |= VC_QCALL; break; + case 'L': vc.flag |= VC_ADJLD; 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; @@ -249,6 +253,7 @@ int bcfview(int argc, char *argv[]) fprintf(stderr, " -H perform Hardy-Weinberg test (slower)\n"); fprintf(stderr, " -N skip sites where REF is not A/C/G/T\n"); fprintf(stderr, " -Q output the QCALL likelihood format\n"); + fprintf(stderr, " -L calculate LD for adjacent sites\n"); fprintf(stderr, " -1 INT number of group-1 samples [0]\n"); fprintf(stderr, " -l FILE list of sites to output [all sites]\n"); fprintf(stderr, " -t FLOAT scaled mutation rate [%.4lg]\n", vc.theta); @@ -259,6 +264,7 @@ int bcfview(int argc, char *argv[]) } b = calloc(1, sizeof(bcf1_t)); + blast = calloc(1, sizeof(bcf1_t)); strcpy(moder, "r"); if (!(vc.flag & VC_VCFIN)) strcat(moder, "b"); strcpy(modew, "w"); @@ -336,6 +342,18 @@ int bcfview(int argc, char *argv[]) if (pr.p_ref >= vc.pref && (vc.flag & VC_VARONLY)) continue; update_bcf1(h->n_smpl, b, p1, &pr, vc.pref, vc.flag); } + if (vc.flag & VC_ADJLD) { // compute LD + double f[4], r2; + if ((r2 = bcf_ld_freq(blast, b, f)) >= 0) { + kstring_t s; + s.m = s.l = 0; s.s = 0; + if (*b->info) kputc(';', &s); + ksprintf(&s, "NEIR=%.3lf", r2); + bcf_append_info(b, s.s, s.l); + free(s.s); + } + bcf_cpy(blast, b); + } if (vc.flag & VC_NO_GENO) { // do not output GENO fields b->n_gi = 0; b->fmt[0] = '\0'; @@ -345,11 +363,10 @@ int bcfview(int argc, char *argv[]) if (vc.prior_file) free(vc.prior_file); if (vc.flag & VC_CALL) bcf_p1_dump_afs(p1); bcf_hdr_destroy(h); - bcf_destroy(b); + bcf_destroy(b); bcf_destroy(blast); vcf_close(bp); vcf_close(bout); if (hash) kh_destroy(set64, hash); if (vc.fn_list) free(vc.fn_list); if (p1) bcf_p1_destroy(p1); return 0; } - diff --git a/bcftools/prob1.c b/bcftools/prob1.c index dc4c16c..e3b0f73 100644 --- a/bcftools/prob1.c +++ b/bcftools/prob1.c @@ -155,8 +155,6 @@ void bcf_p1_destroy(bcf_p1aux_t *ma) } } -#define char2int(s) (((int)s[0])<<8|s[1]) - static int cal_pdg(const bcf1_t *b, bcf_p1aux_t *ma) { int i, j, k; @@ -352,7 +350,7 @@ int bcf_p1_cal(bcf1_t *b, bcf_p1aux_t *ma, bcf_p1rst_t *rst) long double sum = 0.; // set PL and PL_len for (i = 0; i < b->n_gi; ++i) { - if (b->gi[i].fmt == char2int("PL")) { + if (b->gi[i].fmt == bcf_str2int("PL", 2)) { ma->PL = (uint8_t*)b->gi[i].data; ma->PL_len = b->gi[i].len; break; diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl index 0506c5d..eed8766 100755 --- a/bcftools/vcfutils.pl +++ b/bcftools/vcfutils.pl @@ -322,8 +322,8 @@ sub filter4vcf { die(qq/ Usage: vcfutils.pl filter4vcf [options] -Options: -d INT min total depth (given DP or DP4) [$opts{D}] - -D INT max total depth [$opts{d}] +Options: -d INT min total depth (given DP or DP4) [$opts{d}] + -D INT max total depth [$opts{D}] -q INT min SNP quality [$opts{q}] -Q INT min RMS mapQ (given MQ) [$opts{Q}] -1 FLOAT min P-value for strand bias (given PV4) [$opts{1}]