bcftools/main.c

   1 #include <string.h>
   2 #include <stdlib.h>
   3 #include <sys/stat.h>
   4 #include "bcf.h"
   5
   6 int bcfview(int argc, char *argv[]);
   7 int bcf_main_index(int argc, char *argv[]);
   8
   9 #define BUF_SIZE 0x10000
  10
  11 int bcf_cat(int n, char * const *fn)
  12 {
  13         int i;
  14         bcf_t *out;
  15         uint8_t *buf;
  16         buf = malloc(BUF_SIZE);
  17         out = bcf_open("-", "w");
  18         for (i = 0; i < n; ++i) {
  19                 bcf_t *in;
  20                 bcf_hdr_t *h;
  21                 off_t end;
  22                 struct stat s;
  23                 in = bcf_open(fn[i], "r");
  24                 h = bcf_hdr_read(in);
  25                 if (i == 0) bcf_hdr_write(out, h);
  26                 bcf_hdr_destroy(h);
  27 #ifdef _USE_KNETFILE
  28                 fstat(knet_fileno(in->fp->x.fpr), &s);
  29                 end = s.st_size - 28;
  30                 while (knet_tell(in->fp->x.fpr) < end) {
  31                         int size = knet_tell(in->fp->x.fpr) + BUF_SIZE < end? BUF_SIZE : end - knet_tell(in->fp->x.fpr);
  32                         knet_read(in->fp->x.fpr, buf, size);
  33                         fwrite(buf, 1, size, out->fp->x.fpw);
  34                 }
  35 #else
  36                 abort(); // FIXME: not implemented
  37 #endif
  38                 bcf_close(in);
  39         }
  40         bcf_close(out);
  41         free(buf);
  42         return 0;
  43 }
  44
  45 int bcf_main_pwld(int argc, char *argv[])
  46 {
  47         extern double bcf_pair_freq(const bcf1_t *b0, const bcf1_t *b1, double f[4]);
  48         bcf_t *fp;
  49         bcf_hdr_t *h;
  50         bcf1_t **b, *b0;
  51         int i, j, m, n;
  52         double f[4];
  53         if (argc == 1) {
  54                 fprintf(stderr, "Usage: bcftools pwld <in.bcf>\n");
  55                 return 1;
  56         }
  57         fp = bcf_open(argv[1], "rb");
  58         h = bcf_hdr_read(fp);
  59         // read the entire BCF
  60         m = n = 0; b = 0;
  61         b0 = calloc(1, sizeof(bcf1_t));
  62         while (bcf_read(fp, h, b0) >= 0) {
  63                 if (m == n) {
  64                         m = m? m<<1 : 16;
  65                         b = realloc(b, sizeof(void*) * m);
  66                 }
  67                 b[n] = calloc(1, sizeof(bcf1_t));
  68                 bcf_cpy(b[n++], b0);
  69         }
  70         bcf_destroy(b0);
  71         // compute pair-wise r^2
  72         printf("%d\n", n); // the number of loci
  73         for (i = 0; i < n; ++i) {
  74                 printf("%s:%d", h->ns[b[i]->tid], b[i]->pos + 1);
  75                 for (j = 0; j < i; ++j) {
  76                         double r = bcf_pair_freq(b[i], b[j], f);
  77                         printf("\t%.3f", r*r);
  78                 }
  79                 printf("\t1.000\n");
  80         }
  81         // free
  82         for (i = 0; i < n; ++i) bcf_destroy(b[i]);
  83         free(b);
  84         bcf_hdr_destroy(h);
  85         bcf_close(fp);
  86         return 0;
  87 }
  88
  89 int main(int argc, char *argv[])
  90 {
  91         if (argc == 1) {
  92                 fprintf(stderr, "\n");
  93                 fprintf(stderr, "Usage:   bcftools <command> <arguments>\n\n");
  94                 fprintf(stderr, "Command: view      print, extract, convert and call SNPs from BCF\n");
  95                 fprintf(stderr, "         index     index BCF\n");
  96                 fprintf(stderr, "         cat       concatenate BCFs\n");
  97                 fprintf(stderr, "         ld        compute all-pair r^2\n");
  98                 fprintf(stderr, "\n");
  99                 return 1;
 100         }
 101         if (strcmp(argv[1], "view") == 0) return bcfview(argc-1, argv+1);
 102         else if (strcmp(argv[1], "index") == 0) return bcf_main_index(argc-1, argv+1);
 103         else if (strcmp(argv[1], "ld") == 0) return bcf_main_pwld(argc-1, argv+1);
 104         else if (strcmp(argv[1], "cat") == 0) return bcf_cat(argc-2, argv+2); // cat is different ...
 105         else {
 106                 fprintf(stderr, "[main] Unrecognized command.\n");
 107                 return 1;
 108         }
 109         return 0;
 110 }