zal.c

void zalh(int hd1, int hd2, int hd3, int hd4, int ncats)  
{
    int a1, a2, a3 ;  
    ZALLOC(htable, hd1, int ****) ;  
    for (a1=0; a1<hd1; ++a1) { 
     ZALLOC(htable[a1], hd2, int ***) ;  
     for (a2=0; a2<hd2; ++a2)  { 
      ZALLOC(htable[a1][a2], hd3, int **) ;  
      for (a3=0; a3<hd3; ++a3)  { 
        htable[a1][a2][a3] = initarray_2Dint(hd4, ncats, -1) ;
     } 
    }
  }
}

void zalhu(int hd1, int hd2, int hd3, int hd4, int ncats)  
{
    int a1, a2, a3 ;  
    ZALLOC(htableu, hd1, int ****) ;  
    for (a1=0; a1<hd1; ++a1) { 
     ZALLOC(htableu[a1], hd2, int ***) ;  
     for (a2=0; a2<hd2; ++a2)  { 
      ZALLOC(htableu[a1][a2], hd3, int **) ;  
      for (a3=0; a3<hd3; ++a3)  { 
        htableu[a1][a2][a3] = initarray_2Dint(hd4, ncats, -1) ;
     } 
    }
  }
}



void clearhists(HISTS *hpt) 
{
 hpt -> mq = hpt -> bq = hpt -> strand = -1 ; 
 hpt -> pp = -9999 ;   
 hpt -> ignore = NO ;
 hpt -> hitcount = 0 ;
 ivzero(hpt -> hist, 0) ;
}

  
void long2hists(HISTS *hpt, long *a, int ncols) 
{
  double ww[6] ;  
  int kcat ; 
  double y1, y2, z1, z2 ;
  int r , c, t, k ; 
  long *aa ; 
  static double ap, bp ;   

  static int ncall = 0 ; 

  ++ncall ;
  if (ncall==1) { 
   bpars(&ap, &bp, 0.01, 0.01*0.01) ; 
   printf("beta prior for deamination prob: %9.3f %9.3f\n", ap, bp) ;
  }

  hpt -> mq = a[0] ; 
  hpt -> bq = a[1] ; 
  hpt -> strand = a[2] ; 

  t = hpt -> pp = a[3] ; 

  if ((t>readedge) && (t != 99)) { 
   printf("bug: %d\n", ncols) ;  
   printlmat(a, 1, ncols) ;
   fatalx("long2hist bug\n") ;  
  }


  aa = a + 4 ;
  for (k=0; k<16; ++k) { 
   hpt -> hist[k] = aa[k] ;
  }

  kcat = 0 ;  // C T    
  r = 1 ; c = 3 ; 
  y1 = hpt -> hist[4*r+c] + ap ; 
  y2 = hpt -> hist[4*r+r] + bp ;
  z1 = y1/(y1+y2) ;

  kcat = 1 ;  // G A    
  r = 2 ; c = 0 ; 
  y1 = hpt -> hist[4*r+c] + ap ; 
  y2 = hpt -> hist[4*r+r] + bp ;
  z2  = y1/(y1+y2) ;

  if (MIN(z1, z2) < 1.0e-8) fatalx("bad deamin\n") ;  
}

void hists2long(HISTS *hpt, long *a) 
{
  double ww[6] ; 
  a[0] = hpt -> mq  ;
  a[1] = hpt -> bq  ;
  a[2] = hpt -> strand  ;
  a[3] = hpt -> pp  ;

}


//  should probably be in file histinfo.c 
int  getadaptinfo(HISTINFO ***phist, char *iname, int ultmode)  
{

  char ss[MAXSTR] ; 
  int n ; 
  int nrows, ncols ;    
  int  a, b, j, k, x, aa, bb, bx, kk, s, t, pp, pp99, px, jj, tnum ; 
  int a1, a2, a3, a4 ;  
  int  b0, b1, b2, b3, b4, uu[4] ;  
  int maxx, len, r, c, bqlen  ; 
  double y, y1, y2, y3, y4 ; 
  int mq, bq, kill[4] ; 
  long j2[2][16],lx, ls[4], ltot, lll ;
  double xs[4], ydiag ;
  double ynum, yden ;
  double *bcqprob ; 
  double time1, time2 ; 
  char *bamtemp ;
  HISTS histwork, *hwpt ;
  int numdata = 0 ;
  int histempty ;
  char *qsampname ;


  long **data ; 
  double **xx ;  
 
  HISTINFO *hipt ;
  static HISTINFO **histinfo ; 


  int nsplit, maxff=30 ; 
  char *spt[10] ;
  
  long *ww, *w1, *w2  ;
  HISTS *histpt, *h1pt, *h2pt, *hpt, *hpx ; 
  int ***** htab ; 

  hwpt = &histwork ;

  printmem("getadapt") ;
  if (iname == NULL) return -1 ;
  time1 = cputime(1) ;
  n = getfline(ss, iname, MAXSTR) ; 

  nrows = numlines(iname) ;
  ncols = numcolsq(iname)  ;

  // ncols <= 0 means no non-comment (histogram data) line was found: histmake
  // wrote only header/comment lines because zero reads survived filtering for
  // this sample. Treat this as a valid, all-missing histogram rather than a
  // fatal error: nrows/data stay empty, so every later per-read histogram
  // lookup falls through to the existing "missing histogram" (kode -23) path.
  histempty = (ncols <= 0) ;

  if (histempty) {
   printf("## %s :: %s\n", iname, ss) ;
   printf("histfilename: %s has no histogram data (0 reads passed filtering); treating as empty\n", iname) ;
   nrows = 0 ;
   ncols = 0 ;
   data = NULL ;
  }
  else {
   xx = initarray_2Ddouble(ncols, nrows, 0) ;
   nrows = getxxq(xx, nrows, ncols, iname) ;
   if (nrows<0) fatalx("bad format in %s\n", iname) ;

   data = initarray_2Dlong(nrows, ncols, 0) ;
   for (a=0; a<nrows; ++a) {
    for (b=0; b < ncols; ++b) {
     data[a][b] = lrint(xx[b][a]) ;
    }
   }
   printf("## %s :: %s\n", iname, ss) ;
   free2D(&xx, ncols) ;
  }

  maxmq = 50, bq = 50, nump = 10 ; // nump is number of special cols at start and end of read
  qsampname = strdup("???") ;  
  bamtemp  = strdup("???") ; 
  numrgcats = 1 ;

  nsplit = splitup(ss, spt, maxff) ;  
  t = strcmp(spt[0], "###globals") ;
  if (t==0) {
   maxmq = atoi(spt[1]) ;
   maxbq = atoi(spt[2]) ;
   nump = atoi(spt[4]) ; 
   qsampname = strdup(spt[6]) ;
   bamtemp = strdup(spt[8]) ;
   if (ncols>20) numrgcats = atoi(spt[nsplit-1]) ;
  }
  else { 
   printf ("header missing! (defaults applied)") ;
  } 

  if (bamname == NULL) bamname = strdup(bamtemp) ;
  if (sampname != NULL) {
   t = strcmp(qsampname, sampname) ;  
   if (t!=0) printf("*** warning *** sampname does not match\n")  ; 
  }

  printf("numrgcats (getadapt): %d\n", numrgcats) ;  
  fflush(stdout) ; 
  maxx = 2*nump ; 
  hd1 = maxmq + 1 ; 
  hd2 = maxbq + 1 ;
  hd3 = 2 ;  
  hd4 = maxx+1  ;
  if (ultmode == 0) {
   printmem("zalh") ;
   zalh(hd1, hd2, hd3, hd4, numrgcats) ;  
   printmem("zalhx") ;
   htab = htable ;
  }
  if (ultmode == 1) {
   zalhu(hd1, hd2, hd3, hd4, numrgcats) ;  
   htab = htableu ;
  }

  t = 0 ; 
  for (a1=0; a1<hd1; ++a1) {
   for (a2=0; a2<hd2; ++a2) {
    for (a3=0; a3<hd3; ++a3) {
      for (a4=0; a4<hd4; ++a4) {
       ivclear(htab[a1][a2][a3][a4], t, numrgcats) ;  
       ++t ;
      
  }}}} 


  printf("hd dimensions %d %d %d %d\n", hd1, hd2, hd3, hd4) ;
  fflush(stdout) ;

   ZALLOC(histinfo, numrgcats, HISTINFO *) ;

   for (k=0; k<numrgcats; ++k) {
    ZALLOC(histinfo[k], 1, HISTINFO) ;
    hipt = histinfo[k] ; 
    hipt -> maxmq = hipt -> maxbq = -1 ;
    hipt -> nump = nump ;
    hipt -> sampname = NULL ;
    if (sampname != NULL) hipt -> sampname = strdup(sampname) ; 
    hipt -> bamname = strdup(bamname) ; 
    hipt -> bighistlen = 0 ;
    hipt -> rgcat = k ;
    hipt -> isfix = ultmode ;
    strcpy(hipt -> deaminbases, "CGCG") ; // same literal the per-row loop below writes unconditionally; set here too so it's initialized even with zero data rows
   }

    ZALLOC(bglen, numrgcats, int) ;  
    ivclear(bglen, hd1*hd2*hd3*hd4, numrgcats) ;


    bighistlen = nrows ;
     
    if (ultmode == 0) {
     ZALLOC(bighistlist, numrgcats, HISTS **) ;
     bigh = bighistlist ;
    } 
    if (ultmode == 1) {
     ZALLOC(bighistlistu, numrgcats, HISTS **) ;
     bigh = bighistlistu ;
    } 

    for (k=0; k<nrows; ++k) { 
     j = rgnum(data[k], ncols) ;
     if (j>=numrgcats) fatalx("rgnum overflow\n") ; 
    }
     
    for (j=0; j<numrgcats; ++j) { 
     ZALLOC(bigh[j], bglen[j], HISTS *) ;
     bighist = bigh[j] ; 
     for (k = 0; k< bglen[j]; ++k) { 
      ZALLOC(bighist[k], 1, HISTS) ; 
      histpt = bighist[k] ; 
      clearhists(histpt) ;
      hipt = histinfo[j] ; 
      hipt -> bighistlen = bglen[j] ; 
     }
     sprintf(ss, "hists %d loaded\n", j) ;
     printmem(ss) ;
    } 


    maxmq = maxbq = 0 ; 

    ltot = 0 ; 
    
    for (k=0; k<nrows; ++k) { 
     j = rgnum(data[k], ncols) ;
     
     hipt = histinfo[j] ; 
     strcpy(hipt -> deaminbases, "CGCG") ; 
     bighist = bigh[j] ;   
     long2hists(hwpt, data[k], ncols) ;  // 
     x = intsum(hwpt -> hist, 16) ;
     if (x != 0) ++ numdata ;
     a1 = hwpt -> mq ; 
     a2 = hwpt -> bq ; 
     hipt -> maxmq = MAX(hipt -> maxmq, a1) ;
     hipt -> maxbq = MAX(hipt -> maxbq, a2) ;
     if (a1>=hd1) fatalx("overflow 1 %d %d\n", a1, hd1) ;
     if (a2>=hd2) fatalx("overflow 2 %d %d\n", a2, hd2) ;
     maxmq = MAX(maxmq, a1) ;
     maxbq = MAX(maxbq, a2) ;
     a3 = hwpt -> strand ; 
     a4 = hwpt -> pp ; 
     x =  ptoa(a4, nump) ; 
     t = htab[a1][a2][a3][x][j] ; 
     histpt = bighist[t] ;  
     long2hists(histpt, data[k], ncols) ;  // 
   }

  printmem("htabx") ;
  
/**
   bqlen = 4*(maxbq+1) ; 
   ZALLOC(bcqprob, bqlen, double) ;      

   for (j=0; j<numrgcats; ++j) { 
     hipt = histinfo[j] ; 
     bighist = bigh[j] ;   
     len = hipt -> bighistlen ; 
     for (a1=0; a1 <= maxmq; ++a1) { 
      for (a3 =0; a3<2; ++a3) { 
       for (x=0; x<=maxx; ++x) { 
        for (r=0; r<4; ++r) {
        vzero(bcqprob, bqlen) ;
        tnum = 0 ;
        for (t=0; t<len; ++t) {
         histpt = bighist[t] ; 
         if (a1 != histpt -> mq) continue ; 
         if (a3 != histpt -> strand) continue ; 
         a4 = histpt -> pp ; 
         if (x != ptoa(a4, nump)) continue ; 
          bq  = histpt -> bq ; 
          for (c=0; c<4; ++c) { 
            bcqprob[4*bq+c] += (double) histpt -> hist[4*r+c] ;  // r is fixed here 
          }
          ++tnum ; 
        }
 


    }}}}


   }

*/

   free(bglen) ;
   free2Dlong(&data, nrows) ;

   *phist = histinfo ;

  time2 = cputime(1) ; 
  printf("getadaptinfo: %s nonzero rows: %d\n", iname, numdata) ;
  // numdata == 0 with histempty is the expected zero-reads case (handled above).
  // numdata == 0 with data rows present indicates a genuinely malformed histfile.
  if ((numdata == 0) && (!histempty)) fatalx("no data!\n") ;
  printf("getadaptinfo time: %9.3f\n", time2-time1) ;

  printmem("getadaptx") ; 

  return numrgcats ; 
}