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