released/pIRS.old/bwasam/dupsamstat.c

   1 #include <stdlib.h> //EXIT_FAILURE
   2 #include <stdio.h>
   3 //#include <float.h>
   4 #include <math.h>
   5 #include <search.h> //hsearch
   6 #include "bam/sam.h"
   7 #include "bitarray.h"
   8
   9 /*
  10 typedef struct {
  11         int beg, end;
  12         samfile_t *in;
  13 } tmpstruct_t;
  14
  15 // callback for bam_fetch()
  16 static int fetch_func(const bam1_t *b, void *data)
  17 {
  18         bam_plbuf_t *buf = (bam_plbuf_t*)data;
  19         bam_plbuf_push(b, buf);
  20         return 0;
  21 }
  22 // callback for bam_plbuf_init()
  23 static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data)
  24 {
  25         tmpstruct_t *tmp = (tmpstruct_t*)data;
  26         if ((int)pos >= tmp->beg && (int)pos < tmp->end)
  27                 printf("%s\t%d\t%d\n", tmp->in->header->target_name[tid], pos + 1, n);
  28         return 0;
  29 }
  30 */
  31 typedef struct {
  32     void *pBitArray;
  33     uint32_t ChrLen;
  34 } hData_t;
  35
  36 int32_t ChrNum;
  37
  38 int main (int argc, char *argv[]) {
  39         if (argc!=3+1) {
  40                 fprintf(stderr,"Usage: %s <in.sam> <snp.lst> <outprefix>\n",argv[0]);
  41                 exit(EXIT_FAILURE);
  42         }
  43         samfile_t *samIn = samopen(argv[1], "r", NULL);
  44         if (samIn == 0) {
  45                 fprintf(stderr, "Fail to open SAM file %s\n", argv[1]);
  46                 exit(EXIT_FAILURE);
  47         }
  48         bam_header_t *samHeader=samIn->header;
  49         ChrNum = samHeader->n_targets;
  50         hcreate(2*ChrNum);
  51         ENTRY hItem, *phItem;
  52         hData_t *hData = malloc(ChrNum * sizeof(hData_t));
  53     if (hData == NULL) {
  54        fprintf(stderr, "[x]malloc failed\n");
  55        exit(EXIT_FAILURE);
  56     }
  57         for (int32_t i=0;i<ChrNum;++i) {
  58             printf("[%0*d] %21s -> %9u bp\n",(int)(1+log10(ChrNum)),1+i,samHeader->target_name[i],samHeader->target_len[i]);
  59             hItem.key  = samHeader->target_name[i];
  60             hData[i].ChrLen = samHeader->target_len[i];
  61             hData[i].pBitArray = calloc(1, BITNSLOTS(samHeader->target_len[i]) );
  62             hItem.data = &hData[i];
  63         phItem = hsearch(hItem, ENTER);
  64         /* there should be no failures */
  65         if (phItem == NULL) {
  66            fprintf(stderr, "[x]hash entry failed\n");
  67            exit(EXIT_FAILURE);
  68         }
  69
  70         }
  71 /*
  72         if (argc == 2) { // if a region is not specified
  73                 sampileup(samIn, -1, pileup_func, &tmp);
  74         } else {
  75                 int ref;
  76                 bam_index_t *idx;
  77                 bam_plbuf_t *buf;
  78                 idx = bam_index_load(argv[1]); // load BAM index
  79                 if (idx == 0) {
  80                         fprintf(stderr, "BAM indexing file is not available.\n");
  81                         return 1;
  82                 }
  83                 bam_parse_region(samIn->header, argv[2], &ref,
  84                                  &tmp.beg, &tmp.end); // parse the region
  85                 if (ref < 0) {
  86                         fprintf(stderr, "Invalid region %s\n", argv[2]);
  87                         return 1;
  88                 }
  89                 buf = bam_plbuf_init(pileup_func, &tmp); // initialize pileup
  90                 bam_fetch(samIn->x.bam, idx, ref, tmp.beg, tmp.end, buf, fetch_func);
  91                 bam_plbuf_push(0, buf); // finalize pileup
  92                 bam_index_destroy(idx);
  93                 bam_plbuf_destroy(buf);
  94         }
  95 */
  96         samclose(samIn);
  97         for (int32_t i=0;i<ChrNum;++i) {
  98             free(hData[i].pBitArray);
  99         }
 100         free(hData);
 101         return 0;
 102 }