c_cpp/coverages/chrtable.c

   1 #include <stdlib.h>
   2 #include "chrtable.h"
   3 #include <err.h>
   4 #include <zlib.h>
   5
   6 //extern struct ChrData_hash_struct *ChrData;
   7 /*
   8 void add_chr(char *id, uint32_t len) {
   9     struct ChrData_hash_struct *s;
  10
  11     s = malloc(sizeof(struct ChrData_hash_struct));
  12     s->ChrID = malloc(strlen(id)+1);
  13     strcpy(s->ChrID, id);
  14     s->ChrLen = len;
  15     s->Data = malloc(len*sizeof(uint8_t));
  16     HASH_ADD_STR( ChrData, ChrID, s );  // id: name of key field
  17 }
  18 */
  19
  20 char *strlinker(const char * const main, const char * const suffix, char ** outstr) {
  21     char *tmp;
  22     if ((tmp = realloc(*outstr, strlen(main)+strlen(suffix)+1)) == NULL) {
  23         //free(*outstr);
  24         fprintf(stderr, "ERROR: realloc failed\n");
  25     } else {
  26         *outstr = tmp;
  27         strcpy(*outstr,main);
  28         strcat(*outstr,suffix);
  29     }
  30     return *outstr;
  31 }
  32
  33 void inc_depth(int32_t left, int32_t right, uint8_t *ThisDat) {
  34     //if (right <= left) errx(3,"Position Error:[%i,%i].",left,right);
  35     for (int32_t i=left; i<=right; ++i) {
  36         int tmp = ThisDat[i];
  37         if (tmp < UINT8_MAX) {
  38             ++ThisDat[i];
  39         }
  40     }
  41 }
  42
  43 void do_stat(bam1_t *b, const uint16_t overlap, const struct myData *Data) {
  44     uint16_t k = overlap + 1;
  45     //bam1_t *b = balignd;
  46     //bam1_core_t *core = &(b->core);
  47     if (b->core.tid < 0) return;
  48     uint8_t *ThisDat = Data->ChrDat[b->core.tid];
  49     int32_t left = b->core.pos;
  50     int32_t right = bam_calend(&b->core, bam1_cigar(b));
  51     if (left >= right) return;
  52     if (right >= Data->target_len[b->core.tid]) right = Data->target_len[b->core.tid] - 1;
  53     //printf("%i %i %i\n",b->core.tid,left,right);
  54     inc_depth(left, right - k+1, ThisDat);
  55 }
  56
  57 char tmpChr[255];
  58 int do_contig(const uint8_t mindep, const struct myData * Data, const char * const outfile) {
  59     char *tmpoutchar = NULL;
  60     sprintf(tmpChr,"_%u.contig.gz",mindep);
  61     gzFile *fp = gzopen(strlinker(outfile,tmpChr,&tmpoutchar), "wb9");
  62     free(tmpoutchar);
  63 #ifdef DEBUG
  64     sprintf(tmpChr,"_%u.kdepth.gz",mindep);
  65     gzFile *fpd = gzopen(strlinker(outfile,tmpChr,&tmpoutchar), "wb9");
  66     free(tmpoutchar);
  67 #endif
  68     gzprintf(fp,"#minDepth = %u\n",mindep);
  69     int_fast8_t inContig=0;
  70     int32_t lastStart=0;
  71     for (int32_t i=0; i < Data->n_targets; ++i) {
  72         gzprintf(fp,"#-------%s--------\n",Data->target_name[i]);
  73 #ifdef DEBUG
  74         gzprintf(fpd,">%s\n",Data->target_name[i]);
  75 #endif
  76         for (int32_t p=0; p < Data->target_len[i]; ++p) {
  77             if ( (Data->ChrDat[i][p] >= mindep) ) {
  78                 if (!inContig) {
  79                     gzprintf(fp,"%s\t%i\t",Data->target_name[i],p+1);
  80                     inContig=1;
  81                     lastStart = p;
  82                 }
  83             } else {
  84                 if (inContig) {
  85                     gzprintf(fp,"%i\t%i\n",p+1,p-lastStart+1);
  86                     inContig=0;
  87                 }
  88             }
  89 #ifdef DEBUG
  90             gzprintf(fpd,"%i:%u ",p,Data->ChrDat[i][p]);
  91             if (p % 50 == 0) gzprintf(fpd,"\n");
  92 #endif
  93         }
  94         if (inContig) { // if reaching end
  95             gzprintf(fp,"%i\t%i\n",Data->target_len[i],Data->target_len[i] - lastStart);
  96             //inContig=0;
  97         }
  98 #ifdef DEBUG
  99         gzprintf(fpd,"\n");
 100 #endif
 101     }
 102     gzclose(fp);
 103 #ifdef DEBUG
 104     gzclose(fpd);
 105 #endif
 106     return 0;
 107 }