bcftools/vcf.c

   1 #include <zlib.h>
   2 #include <stdlib.h>
   3 #include <stdio.h>
   4 #include <string.h>
   5 #include "bcf.h"
   6 #include "kstring.h"
   7 #include "kseq.h"
   8 KSTREAM_INIT(gzFile, gzread, 4096)
   9
  10 typedef struct {
  11         gzFile fp;
  12         FILE *fpout;
  13         kstream_t *ks;
  14         void *refhash;
  15         kstring_t line;
  16 } vcf_t;
  17
  18 bcf_hdr_t *vcf_hdr_read(bcf_t *bp)
  19 {
  20         kstring_t meta, smpl;
  21         int dret;
  22         vcf_t *v;
  23         bcf_hdr_t *h;
  24         if (!bp->is_vcf) return 0;
  25         h = calloc(1, sizeof(bcf_hdr_t));
  26         v = (vcf_t*)bp->v;
  27         v->line.l = 0;
  28         memset(&meta, 0, sizeof(kstring_t));
  29         memset(&smpl, 0, sizeof(kstring_t));
  30         while (ks_getuntil(v->ks, '\n', &v->line, &dret) >= 0) {
  31                 if (v->line.l < 2) continue;
  32                 if (v->line.s[0] != '#') return 0; // no sample line
  33                 if (v->line.s[0] == '#' && v->line.s[1] == '#') {
  34                         kputsn(v->line.s, v->line.l, &meta); kputc('\n', &meta);
  35                 } else if (v->line.s[0] == '#') {
  36                         int k;
  37                         char *p, *q, *r;
  38                         for (q = v->line.s, p = q + 1, k = 0; *p; ++p) {
  39                                 if (*p == '\t' || *(p+1) == 0) {
  40                                         r = *(p+1) == 0? p+1 : p;
  41                                         if (k >= 9) {
  42                                                 kputsn(q, r - q, &smpl);
  43                                                 kputc('\0', &smpl);
  44                                         }
  45                                         q = p + 1; ++k;
  46                                 }
  47                         }
  48                         break;
  49                 }
  50         }
  51         kputc('\0', &meta);
  52         h->name = 0;
  53         h->sname = smpl.s; h->l_smpl = smpl.l;
  54         h->txt = meta.s; h->l_txt = meta.l;
  55         bcf_hdr_sync(h);
  56         return h;
  57 }
  58
  59 bcf_t *vcf_open(const char *fn, const char *mode)
  60 {
  61         bcf_t *bp;
  62         vcf_t *v;
  63         bp = calloc(1, sizeof(bcf_t));
  64         v = calloc(1, sizeof(vcf_t));
  65         bp->is_vcf = 1;
  66         bp->v = v;
  67         if (strchr(mode, 'r')) {
  68                 v->fp = strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
  69                 v->ks = ks_init(v->fp);
  70         } else if (strchr(mode, 'w'))
  71                 v->fpout = strcmp(fn, "-")? fopen(fn, "w") : fdopen(fileno(stdout), "w");
  72         return bp;
  73 }
  74
  75 void bcf_hdr_clear(bcf_hdr_t *b);
  76
  77 int vcf_close(bcf_t *bp)
  78 {
  79         vcf_t *v;
  80         if (bp == 0) return -1;
  81         if (bp->v == 0) return -1;
  82         v = (vcf_t*)bp->v;
  83         if (v->fp) {
  84                 ks_destroy(v->ks);
  85                 gzclose(v->fp);
  86         }
  87         if (v->fpout) fclose(v->fpout);
  88         free(v->line.s);
  89         free(v);
  90         free(bp);
  91         return 0;
  92 }
  93
  94 int vcf_hdr_write(bcf_t *bp, const bcf_hdr_t *h)
  95 {
  96         vcf_t *v = (vcf_t*)bp->v;
  97         int i;
  98         if (v == 0 || v->fpout == 0) return -1;
  99         fwrite(h->txt, 1, h->l_txt, v->fpout);
 100         fprintf(v->fpout, "#CHROM\tPOS\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT");
 101         for (i = 0; i < h->n_smpl; ++i)
 102                 fprintf(v->fpout, "\t%s", h->sns[i]);
 103         fputc('\n', v->fpout);
 104         return 0;
 105 }
 106
 107 int vcf_read(bcf_t *bp, bcf1_t *b)
 108 {
 109         int dret;
 110         vcf_t *v = (vcf_t*)bp->v;
 111         v->line.l = 0;
 112         if (ks_getuntil(v->ks, '\n', &v->line, &dret) < 0) return -1;
 113         return v->line.l + 1;
 114 }
 115
 116 int vcf_test(int argc, char *argv[])
 117 {
 118         bcf_t *bp, *bpout;
 119         bcf_hdr_t *h;
 120         bp = vcf_open(argv[1], "r");
 121         bpout = vcf_open("-", "w");
 122         h = vcf_hdr_read(bpout);
 123         vcf_hdr_write(bpout, h);
 124         vcf_close(bp);
 125         return 0;
 126 }