* 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 2dbcf594ee84c2fae9c183bdd43e3c1480973b3f..b6ab9726a6efa25e1554f7cb28aab4d1fc4a9600 100644 (file)
--- 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 35c3d943e04a62fa65387394d141946d7e16b47b..3531647e18775cd5049954385a1022d39a17ec6f 100644 (file)
--- 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 72861113d5bac07b66d7014f41a8c95b888e143d..36578956bbae06d4c216c121a2ddb3a94cad94d8 100644 (file)
--- 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 1eacd4f731f1e35866a024d024ffe24cf4d8b49b..aa3e591a3d859f2ee834dcff6715d56920d99ca0 100644 (file)
--- 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 dc4c16ccc7c24edafd39ec0de193eb618435beb8..e3b0f7344d096f22c37048b2131219378bcb7ba5 100644 (file)
--- 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 0506c5dbacb4e4480eebf4eab5990787d423c7f9..eed8766464e9b63e08268a97b146d3464705a042 100755 (executable)
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -322,8 +322,8 @@ sub filter4vcf {
   die(qq/
 Usage:   vcfutils.pl filter4vcf [options] <in.vcf>
 
-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}]