histmake.c

#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <math.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <libgen.h>

#include <nicklib.h>
#include <globals.h>
#include <getpars.h>
#include <nicksam.h>

#include "admutils.h"
#include "mcmcpars.h"
#include "mcio.h"  
#include "popsubs.h" 

#include "cigsubs.h" 

#define WVERSION   "1850" 

// printcount added
// meandepth added
// depth > 100 counted as 100
// rgonly added
// single-stranded UDG added
// badrgname added 
// rgcount added
// printbadrgs added 
// printrghit  added
// rgnone count added
// singlebam mode compulsory (bamname:   sampname:) 
// snp alleles must be set valid
// readdepth added (big output file) ;
// maxfullpops  1000 
// very long reads trapped 
// MAXPOPS 5000
// code cleaned for O2
// new db format (4 fields) upward compatible
// readedge added
// rgcats added
// polymask added  (shop 1000G polymorphic fa)
// rgsplit added
// killcpg  added  
// cpgonly  added  
// minlen, maxlen added
// rgcount added
// printread added  (use with maxread) 
// deaminhistname added (Daniel request)
// chrom 90 (MT) supported
// countbasepos added
// cadjustname added
// maxcl fixed (calcudg) 
// empty hists not written 
// mqclip 
// deaminhistname fixed for indels
// more careful memory allocation  
// processdups added
// countbases added (used by damagescore in drfancy )
// forcereg added
// killsnps added
/** 
 For single chromosome organisms with non-integer chromosome names 
forcereg: <name> 
allows one to pretend the chromosome is 1 for bam and reference.fa processing
It would be easy to add a feature to map chromosomes 1,2,... to names. 
Similar change made to pullit. 

killsnps: YES applies a mask to reference from snpfile.  Note that histmake output
will now depend on snpname: <snpfile>

*/


int ***LC, ***RC ;

long ***cbases ;   // numrg X 2 X 12 

int mapreg = YES ; 

#define MAXFL  50   
#define MAXSSS  8192
#define MAXPOPS 5000
#define MAXRG 50000 

#define MAXL 150 
#define MINL 8 

int pulldownmode = NO ; 
int pullmode = NO ;
int printtable = NO ;
int *rgpops, numrg ;
char **iublist = NULL ; 
int flipstrand = NO ;  // debug/test only
int printrd = NO ;
int posbasestats = NO ; 
char *cadjustname = NULL ;
int cadjstackbase = NO ; 
//  YES => base dist is stack dependent
int ultima = NO ;  
int bqclip = 999 ;  
int mqclip = 999 ;  
long xindels[2] ; 
int pubcigar = 0 ;
long badcigar = 0 ;
int lomap = 10, loqual = 20, ispair = NO  ;

int *threshlo, *threshhi, *quallo ;
char *iubpops = NULL ;
int numiub = 0 ;
int printcount = NO ;
int usechimp = NO ;
int rgsplit = NO ;

char *dstag = "ds:f:" ; 
char *rctag = "rc:i:" ; 

char *iubfile = NULL ;
char *dbfile = NULL ;
char *dbbam = NULL ;
char *samxtra = NULL ;
char *readbam = NULL ; 
char *bamname = NULL, *sampname = NULL ;
char *readdepthname = NULL ;
char *samtoolspath="samtools" ; 
char *deaminhistname = NULL ;
FILE *deaminfile = NULL;  
int  singlebam = YES ;
// int maxread = 1000*1000*1000 ; // debug
long maxread = -1 ;
char xregname[10], yregname[10] ;

char *parname = NULL ;
char *xtag = NULL ;
char *oldlogfilename = NULL ; 


int rgonly = NO ; // yes prints rg group and quits
int usepmdscore = NO ; // on read.  Just effects readbam
int flatbase = NO ;

int processdups = NO ;

double pmdlo = -999 ;
double pmdhi = 999 ;

int minlen = MINL ; 
int maxlen = MAXL ; 

char xregname[10], yregname[10], mtregname[10] ;

char *udgmode = NULL ;
int udgtype = -1 ;
int calcudgtype = YES ;
int udgverbose = NO ;

extern int checksizemode ;
// the variables below should not be needed
char *omode = "eigenstrat" ;     // compulsory as we only support 
extern int packmode ;            //!< flag - input {is not,is} in packed mode
extern char *packgenos ;         //!< packed genotype data (packit.h)
extern char *packepath ;
extern long packlen;             //!< allocated size of packgenos data space
extern long rlen;                //!< number of bytes in packgenos space that each SNP's data occupies

char *trashdir = "/var/tmp" ;
extern int verbose  ;
int qtmode = NO ;
Indiv **indivmarkers;
SNP **snpmarkers = NULL ;
int numsnps, numindivs ; 
int hasrggroup = YES ;
int basewindow = 0 ;

long nreject, nindelx, nindel  ;
int readsumm = 0, readsummc = 0 ;
long cpgmiss = 0 ;

char  *genotypename = NULL ;
char  *snpname = NULL ;

char *indoutfilename = NULL ;
char *snpoutfilename = NULL ;
char  *genooutfilename = NULL ;

char  *indivname = NULL ;

int checkmode = NO ;

char **samplist ; 
int nsamps ;

int minchrom = 1 ;
int maxchrom = 24 ;
int xchrom = -1 ;

char *reffasta = NULL ;

char *outputname = NULL ;
FILE *ofile ;

double fakespacing = 0.0 ;
char  unknowngender = 'U' ;

char *regname = NULL ;
char **iublist ;
char *regstring = NULL ; 
int reglen ;
char *chimpstring  = NULL; 
int  chimplen ;


char *rstring ; 

char *poplistname = NULL ;
char **poplist ;
char **rawpoplist ;
char **fullpoplist ;

char *badrgname = NULL ; 
    
char **bad ; 
int nbadrg = 0 ;

int npops=0, nfpops = 0, ret ; 
int seed = 0  ; 

char *fasta ;

char *histfilename = NULL ;
char *ref[1] ;  // ref seq name
char *chimpref[1] ;  // chimp seq name

char *polyref[1] ;  // poly seq name
char *polystring  = NULL; 
char *polyname = "polymask" ;
int  polylen ;
int polymask = YES ;  
int killcpg = NO ;
int cpgonly = NO ;
int stopcount  = NO ;

READ  **readlist, *readpt, *readptx, *oldreadpt ; 
THRESH *pthresh, ptx ;

  HIST ***bighistlist ;  
  HIST ***bighistlistu ;  
  HIST **bighist, *histpt ; 
  int maxbighistlen = 500000, bighistlen ;  
  int hd1 = 101, hd2=101, hd3 =2, hd4 = 50 ;
  int ****htable ; 
  int *ptoax, ppwk[200], numpp, readedge = 10  ; 
  int maxmq=0, maxbq = 0 ;

char *rgcatsname = NULL ;
int numrgcats = 1, numrgc = 0 ;  
int rgnone = 0 ;
char **rgclist ; 
int *rgcats, *rgcount, rgmiss = 0  ; 

char *refname = "Href" ;
int *cat2rg, *rg2cat ; 
char *forcereg = NULL ; 
int killsnps = NO ;

void readcommands(int argc, char **argv) ;
long setgenoblank (SNP **snpmarkers, int numsnps, int numindivs)   ;
int readfa1(char *faname, char **pfasta, int *flen) ;
char *myfai_fetch(faidx_t *fai, char *reg, int  *plen) ;
int  mksamplist(char **samplist, Indiv **indivmarkers, int numindivs) ;
void loadreg(char *reg) ;   
void mkpop (char **rawpops, char **pops, int n) ;
void calccount(int ***countarr, int firstsnp, int lastsnp) ;
void fixreg(int chrom) ;
int getgt(int *aa) ;
int stripindivs(Indiv **indivmarkers, int numindivs) ;

int getfalist(char **poplist, int npops, char *dbfile, char **iublist)  ;
int checkread(READ *readpt)  ;
int setudgtype(char *mode) ; 
void clearhist(HIST *hpt)  ;
void counthist(READ *readpt, HIST **bighist)  ;
void clear4D(int ****table, int x1, int x2, int x3, int x4, int val) ;  
void dumphist (char *histfilename, HIST *** bigh, int ultmode) ;
void zalh(int hd1, int hd2, int hd3, int hd4)  ;
void makergcats(char *rrr, char *bamname) ;
int qcpg1(int pos, char *regstring, int reglen) ;
int countsam(int first, int last) ; 
int callsam() ; 
int ptox(int p)  ;
void crackx(long *bb, HIST *hpt)   ;
char *udgm(int ds, int tt) ;
void calcudg(HIST **bighist, int *stype, int *utype) ;
void deaminhist(HIST **bighist, FILE *fff) ;
void countbasepos () ;
void  setthreshx(THRESH *pthresh)  ;
int testconchi(double *a, int m, int n) ;
 
FILE *samstream ;

int main(int argc, char **argv)
{

  int i, j, k, d, n, kk, g, s, t, tt, ccheck ; 
  int pos, isnp, jsnp, gval ;
  SNP *cupt ;
  Indiv *indx ;
  char *chimpname = "Chimp" ;

  int numvind, nignore, numrisks = 1 ;
  int chrom, lastchrom = -1 ;  
  int firstsnp, rlen, mlen  ;
  int lastsnp, lensnp ; 
  int ***countarr, depth ;
  char ss[MAXSTR] ;
  int *dsum, *dnum, *dcover ;
  int isok, t1, t2 ; 
  double meandepth, ymem ;
  FILE *readdfff = NULL, *fjunk ;

  int mq, bq, pp, ww[100] ; 
  int a, b ;
  int d1, d2, d3, d4 ;  
  int a1, a2, a3, a4 ;  
  double **xx ;  
  int **snppos, *fpos, *lpos ; 
  double y0, y1, yy ;

  cputime(0) ;
  calcmem(0) ;

  ofile = stdout; 
  packmode = YES ;
  ref[0] = polyref[0] = NULL ; 
  readcommands(argc, argv) ;

  if (polyref[0] != NULL) polymask = YES ;

   if (parname == NULL) { 
    printf("version: %s\n", WVERSION) ; 
    return 0 ;
   }


  printf("executable: %s\n", argv[0]) ; 
  printf("version: %s\n", WVERSION) ; 
  printnl() ;
  printcmdline(argc, argv)  ;
  printf("cwd: %s\n", getcwd(ss, MAXSTR)) ;
  printnl() ;
  setdump(coredump) ;

  setref(&reffasta, iubfile, "Href") ;     
  ref[0] = strdup(reffasta) ;
  printf("reference: %s\n", reffasta) ; 
 
  if (histfilename == NULL) fatalx("no histfile!\n") ;
  printf("histfilename: %s\n", histfilename) ;
  if (forcereg) printf("forcereg: %s\n", forcereg) ;
  if (killsnps) {
   numsnps = getsnps(snpname, &snpmarkers, 0,  NULL, &nignore, numrisks) ;
  } 

  hd1 = MIN(mqclip+1, hd1) ;
  hd2 = MIN(bqclip+1, hd2) ;
  hd4 = 2*readedge + 1 ;

 maxbighistlen = hd1*hd2*hd3*hd4 ;
 ++maxbighistlen ;

  if (bqclip < 900) printf("bqclip: %d\n", bqclip) ;
  if (mqclip < 900) printf("mqclip: %d\n", mqclip) ;
  if (udgverbose) printf("udgverbose set\n") ;

  nreject = nindelx = nindel = 0 ;
  if (cadjustname != NULL) posbasestats = YES ;

  if (basewindow > 0) { 
   t = basewindow % 2 ; 
   if (t==0) ++basewindow ;  // force odd
  }

  if (samxtra != NULL) printf("samxtra: %s\n", samxtra) ; 

  readdepthname = NULL ;

  printf("histmake: version: %s\n", WVERSION) ; 
  if (cpgonly) killcpg = NO ;
  if (outputname != NULL) openit(outputname, &ofile, "w") ;

  singlebam = YES ;
  if (ultima) printf("ultima mode set!\n") ;

  if (seed == 0) seed = filehash(parname) ; 
  SRAND(seed) ; 
  printf("seed: %d\n", seed) ;

  if (reffasta != NULL) printf("reference fasta: %s\n", reffasta) ;

  if (numchrom<0) numchrom = 22 ;  
  if (numchrom != 22)  printf("numchrom: %d\n", numchrom) ;
  sprintf(xregname, "%d", numchrom+1) ;  
  sprintf(yregname, "%d", numchrom+2) ;  

  if (readdepthname != NULL)  openit(readdepthname, &readdfff, "w") ;
  minlen = MAX(minlen, 2*readedge) ; 

  numindivs = 1 ; 

  xindels[0] = xindels[1] = 0 ;
  
  ZALLOC(samplist, numindivs, char *) ;
  if (sampname == NULL) fatalx("no sampname!\n") ;
  samplist[0] = strdup(sampname) ;

  nsamps = 1 ;
  poplist = samplist ;
  if ((singlebam) && (nsamps != 1)) fatalx("singlebam set but nsamps !=  1\n") ;
  if ((singlebam) && (sampname == NULL)) fatalx("singlebam set but no sampname\n") ; 
  if ((bamname == NULL) && (sampname != NULL)) fatalx("sampname set but not bamname\n") ; 

  t = openit_trap(bamname, &fjunk, "r") ; 
  if (t==NO) fatalx("failed to open %s\n", bamname) ;
  fclose(fjunk) ;


   sprintf(xregname, "%d", numchrom+1) ;  
   sprintf(yregname, "%d", numchrom+2) ;  
   sprintf(mtregname, "%d", 90) ;  
  
  
  
  printf("samplist:\n") ;
  printstrings(samplist, nsamps) ;
  
  
  printf("numsnps: %d  numindivs:  %d\n", numsnps, numindivs) ;
  fflush(stdout) ;

   if (polymask) {
    getfalist(&polyname, 1, iubfile, polyref) ;  // ref[0] is now reference fasta ;  
    printf("polymask: %s\n", polyref[0]) ; 
   }
  


 ZALLOC(LC, 2, int **) ;
 ZALLOC(RC, 2, int **) ;
 for (s=0; s <= 1; ++s) { 
  LC[s] = initarray_2Dint(4,4,0) ;
  RC[s] = initarray_2Dint(4,4,0) ;
 }


   npops = 1 ;
   ZALLOC(poplist, npops, char *) ;
   ZALLOC(rawpoplist, npops, char *) ;
   ZALLOC(fullpoplist, MAXPOPS, char *) ;

   if (sampname == NULL) {
    fatalx("no sampname!\n") ;
    npops = loadlist(rawpoplist, poplistname) ;
   }
   else {
    fullpoplist[0] = poplist[0] = rawpoplist[0] = strdup(sampname) ; 
   }

   ZALLOC(dsum, numindivs, int) ;
   ZALLOC(dnum, numindivs, int) ;
   ZALLOC(dcover, numindivs, int) ;

// must fix indivmarkers for missings.
  pthresh = &ptx ;

  setthreshx(pthresh) ;    

  ZALLOC(rgpops, MAXRG, int) ;

  nbadrg = 0 ;
  if (badrgname != NULL) { 
   nbadrg = numlines(badrgname) ; 
   ZALLOC(bad, nbadrg, char *) ; 
   nbadrg = loadlist(bad, badrgname) ;
  }

   numrg = 0 ;

  if (rgsplit) {  
   rgcatsname = mytemp("rggroups") ; 
   if (bamname == NULL) fatalx("no bamname!\n") ;
   makergcats(rgcatsname, bamname) ;
  }
  numrgcats = 1 ;  
  if (rgcatsname != NULL) {  
   t = numlines(rgcatsname) ;  
   ZALLOC(rgclist, t, char *) ; 
   ZALLOC(rgcats, t, int) ;  
   rgmiss = 0 ;
   xx = initarray_2Ddouble(1, t, -1) ; 
   numrgc = getxxnames(&rgclist, xx, t, 1, rgcatsname) ;
   fixit(rgcats, xx[0], numrgc) ; 
    
   free2D(&xx, 1) ; 
   ivmaxmin(rgcats, numrgc, &numrgcats, NULL) ;  
   ++numrgcats ;  
   ZALLOC(rgcount, numrgcats, int) ;  
   ZALLOC(cat2rg, numrgcats, int) ;  
   ZALLOC(rg2cat, numrgc, int) ;  
   printnl() ; 
   printf("rgcats: %d\n", numrgcats) ;  
   hasrggroup = YES ;
   for (k=0; k< numrgc; ++k) { 
    t = rgcats[k] ;
    cat2rg[t] = k ;
    rg2cat[k] = t ;
   }
   for (k=0; k< numrgc; ++k) { 
    printf("%30s %3d\n", rgclist[k], rgcats[k]) ;
   }
   printnl() ; 
  }
  else {
   hasrggroup = NO ;
   ZALLOC(rgcount, numrgcats, int) ;  
   ZALLOC(cat2rg, numrgcats, int) ;  
   ZALLOC(rg2cat, numrgc, int) ;  
  }

  if (rgonly) { 
   printf("rgonly set\n") ;
   printf("##end of histmake\n") ;
   return 0 ;
  }

  countarr = NULL; lensnp = -1 ;
  t = 3 ;

  ZALLOC(readlist, t, READ *) ; 
  ZALLOC(readlist[0], 1, READ) ; 
  readpt = readlist[0] ;
  cleanread(readpt) ;

  ZALLOC(readlist[1], 1, READ) ; 
  readptx = readlist[1] ;
  cleanread(readptx) ;
  
  ZALLOC(readlist[2], 1, READ) ; 
  oldreadpt = readlist[2] ;
  cleanread(oldreadpt) ;
  

    t = 2*readedge + 1 ;  if (t>50) fatalx("read window too big\n") ;

    ZALLOC(bighistlist, numrgcats, HIST **)  ; 
    ZALLOC(cbases, numrgcats, long **)  ; 

    for (j=0; j<numrgcats; ++j) { 

     cbases[j] = initarray_2Dlong(2, 12, 1) ;
     
     ZALLOC(bighistlist[j], maxbighistlen, HIST *) ;
     bighist = bighistlist[j] ; 

      for (k = 0; k< maxbighistlen; ++k) { 
       ZALLOC(bighist[k], 1, HIST) ; 
       histpt = bighist[k] ; 
       clearhist(histpt) ;
      }
    }

    ZALLOC(bighistlistu, numrgcats, HIST **)  ; 
    for (j=0; j<numrgcats; ++j) { 
     
     ZALLOC(bighistlistu[j], maxbighistlen, HIST *) ;
     bighist = bighistlistu[j] ; 

      for (k = 0; k< maxbighistlen; ++k) { 
       ZALLOC(bighist[k], 1, HIST) ; 
       histpt = bighist[k] ; 
       clearhist(histpt) ;
      }
    }
 
    t = 0 ; 
    ivclear(ppwk, -9999, 200) ;
    ptoax = ppwk+100 ; // index can be negative

    for (k=0; k<readedge; ++k) { 
     pp = k + 1 ; 
     ww[t] = pp ;
     ++t ;  
    }
    for (k=0; k<readedge; ++k) { 
     pp = -(k + 1) ; 
     ww[t] = pp ;
     ++t ;  
    }

    pp = 99 ;  
    ww[t] = pp ;
    ++t ;  

    numpp  = t ;

    for (a=0; a<numpp; ++a) { 
     pp = ww[a] ;  
     ptoax[pp] = a ;  
    }

    hd4 = 2*readedge+2 ; 
    zalh(hd1, hd2, hd3, hd4) ;

    t = 0 ;
    clear4D(htable, hd1, hd2, hd3, hd4, -1) ;
    for (mq=0; mq<hd1; ++mq) { 
     for (bq=0; bq<hd2; ++bq) { 
      for (s=0; s<=1; ++s) { 
       for (a=0; a<numpp; ++a) { 
        pp = ww[a] ;  
        htable[mq][bq][s][a] = t ; 
        for (j=0; j<numrgcats; ++j) {
         bighist = bighistlist[j] ; 
         histpt = bighist[t] ; 
         histpt -> mq = mq ; 
         histpt -> bq = bq ; 
         histpt -> strand = s ; 
         histpt -> pp = pp ;
         bighist = bighistlistu[j] ; 
         histpt = bighist[t] ; 
         histpt -> mq = mq ; 
         histpt -> bq = bq ; 
         histpt -> strand = s ; 
         histpt -> pp = pp ;
        } 
        ++t ;
        if (t>=maxbighistlen) fatalx("overflow (bighist)\n") ;    
       }
      }
     }
    }
    bighistlen = t ;
    printf("bighistlen: %6d max: %6d\n", bighistlen, maxbighistlen) ;

   if (xchrom > 0) minchrom = maxchrom = xchrom ; 

   for (chrom=minchrom; chrom <= maxchrom; ++chrom) { 

    sprintf(ss, "%d", chrom) ;
    if (forcereg) strcpy(ss, forcereg) ;
    loadreg(ss) ;
    fixreg(chrom) ; 

    printf("calling sam\n") ; fflush(stdout) ; 
    callsam() ; 
    printf("called sam\n") ; fflush(stdout) ; 
    n = countsam(1, reglen) ; 
    printf("reg: %s  reads: %d\n", regname, n) ; fflush(stdout) ;

    printf("chrom %d processed\n", chrom) ;
    fflush(stdout) ;

  }

 if (numrg>0) {
  printbadrgs() ;
  printrghit() ;
  if (rgnone>0) printf("*** reads with no RG group: %d\n", rgnone) ; 
 }

  countbasepos()  ;

  deaminfile = stdout ; 
  if (deaminhistname != NULL) openit(deaminhistname, &deaminfile, "w") ; 


  for (k=0; k<numrgcats; ++k) {  
   for (j=0; j<=1; ++j) { 
    printf("cbases rcat: %d strand: %d\n", k, j) ; 
    printlmat(cbases[k][j], 3, 4) ; 
   }
  }

 dumphist(histfilename, bighistlist, ultima) ; 

 if (deaminhistname  != NULL) { 
  fclose(deaminfile) ;
  deaminfile = NULL ;
  deaminhistname = NULL ; 
 }


  printf("badcigars: %ld\n", badcigar) ;
  y0 = xindels[0] ; 
  y1 = xindels[1] ; 
  yy = y1 / (y0 + y1 + 1.0e-6) ; 

  printf("indel count (reject, no indel, indel)" ) ;
  printf(" %ld", nreject) ;
  printf(" %ld", nindelx) ;
  printf(" %ld", nindel) ;
  if (cpgonly || killcpg) printf(" cpgmiss: %ld", cpgmiss) ;
  printnl() ;
  printf("indel rate: %9.3f\n", yy) ; 

  ymem = calcmem(1)/1.0e6 ; 
  printf("##end of histmake: %12.3f seconds cpu %12.3f Mbytes in use\n", cputime(1), ymem) ;

  return 0 ;

  ymem = calcmem(1)/1.0e6 ; 
  printf("##end of histmake: %12.3f seconds cpu %12.3f Mbytes in use\n", cputime(1), ymem) ;

  return 0 ;
}
int getgt(int *aa) 
{
  int depth ; 
  double pref ;

  depth = intsum(aa, 2) ;
  if (depth==0) return -1 ;

  if (aa[0]==aa[1]) return 2*ranmod(2) ;


  if (aa[0]==0) return 0 ;
  if (aa[1]==0) return 2 ;

  pref = (double) aa[0] / (double) depth ; 
  if (prob1(pref)==YES) return 2 ; 
  return 0 ;
}

void fixreg(int chrom) 
// reference set OK except on target snps
{
   char *treg = NULL ;
   int k, len, pos ;
   SNP *cupt ; 

// lower case used to mask bases not in snpname
   len = strlen(regstring) ;

   for (k=0; k<numsnps; ++k) { 
    cupt = snpmarkers[k] ; 
    if (cupt -> chrom != chrom) continue ; 
    pos = snpmarkers[k] -> physpos ;
    if (pos <= 0) continue ;
    if (pos >= len) continue ;
    regstring[pos-1] = '?' ;
   }
}
  
void readcommands(int argc, char **argv) 

{
  int i,haploid=0;
  phandle *ph ;
  char str[5000]  ;
  char *tempname ;
  char *tparams ; 
  int n ;
  int targc; 
  char *targv[5] ; 

  while ((i = getopt (argc, argv, "r:p:t:cvVs")) != -1) {

    switch (i)
      {

      case 'p':
	parname = strdup(optarg) ;
	break;

      case 't':
	xtag = strdup(optarg) ;
	break;


      case 'v':
	printf("version: %s\n", WVERSION) ; 
	break; 

      case 'c':
	checkmode = YES ;
	break; 


      case 'V':
	verbose = YES ;
	break; 

      case 's':
	seed = seednum() ;  // random seed;  No exact replay 
	break;

      case 'r':
	oldlogfilename = strdup(optarg) ;
        randomname(&tempname) ;
        tparams = mytemp(tempname) ;
        mkparamfile(oldlogfilename, tparams) ; 
        parname = strdup(tparams) ;
        break ; 

      case '?':
	printf ("Usage: bad params.... \n") ;
	fatalx("bad params\n") ;
      }
  }
         
   if (parname == NULL) { 
    printf("no parameters!\n") ;
    return ;
   }
   printf("##parameter file: %s\n", parname) ;
   ph = openpars(parname) ;
   dostrsub(ph) ;
   fixpars(ph, xtag) ;
   fixpars(ph, basename(argv[0])) ;

   getstring(ph, "bamname:", &bamname) ;  
   getstring(ph, "cramname:", &bamname) ;  
   getstring(ph, "sampname:", &sampname) ;  
   getstring(ph, "refname:", &refname) ;  
   getstring(ph, "referencefile:", &reffasta) ;  
   getint(ph, "minchrom:", &minchrom) ;
   getint(ph, "maxchrom:", &maxchrom) ;
   getint(ph, "chrom:", &xchrom) ;
   getint(ph, "pmdscore:", &usepmdscore) ;
   getint(ph, "rgonly:", &rgonly) ;
   getint(ph, "flipstrand:", &flipstrand) ;
   getlong(ph, "maxread:", &maxread) ;
   getint(ph, "seed:", &seed) ;
   getint(ph, "numchrom:", &numchrom) ;
   getint(ph, "calcudgtype:", &calcudgtype) ;
   getint(ph, "printread:", &printrd) ;
   getint(ph, "printrd:", &printrd) ;
   getint(ph, "posbasestats:", &posbasestats) ;
   getint(ph, "flatbase:", &flatbase) ;
   getint(ph, "ultima:", &ultima) ;
   getint(ph, "indelaware:", &ultima) ;
   getint(ph, "bqclip:", &bqclip) ;
   getint(ph, "mqclip:", &mqclip) ;
   getint(ph, "pisubcigar:", &pubcigar) ;
   getint(ph, "lomap:", &lomap) ;
   getint(ph, "loqual:", &loqual) ;
   getint(ph, "pair:", &ispair) ;
   getint(ph, "ispair:", &ispair) ;
   getint(ph, "readsumm:", &readsumm) ;
   getint(ph, "processdups:", &processdups) ;
   getint(ph, "coredump:", &coredump);
   getint(ph, "cadjstackbase:", &cadjstackbase) ;
   getint(ph, "udgverbose:", &udgverbose) ;

   getdbl(ph, "pmdlo:", &pmdlo) ;
   getdbl(ph, "pmdhi:", &pmdhi) ;

   getstring(ph, "fadatabase:", &iubfile) ;
   getstring(ph, "samxtra:", &samxtra) ;
   getstring(ph, "samtools:", &samtoolspath) ;
   getstring(ph, "samtoolspath:", &samtoolspath) ;
   getstring(ph, "readbam:", &readbam) ;
   getstring(ph, "badrgname:", &badrgname) ;
   getstring(ph, "readdepthname:", &readdepthname) ;
   getstring(ph, "ref:", &reffasta) ;  
   getstring(ph, "polyname:", &polyname) ;
   getstring(ph, "polymaskfile:", &polyref[0]) ;
   getstring(ph, "histfilename:", &histfilename) ;
   getstring(ph, "deaminhistname:", &deaminhistname) ;
   getstring(ph, "deaminhist:", &deaminhistname) ;
   getstring(ph, "rgcats:", &rgcatsname) ;
   getstring(ph, "rgcatsname:", &rgcatsname) ;
   getint(ph, "readedge:", &readedge) ;
   getint(ph, "killcpg:", &killcpg) ;
   getint(ph, "cpgonly:", &cpgonly) ;
   getint(ph, "polymask:", &polymask) ;
   getint(ph, "rgsplit:", &rgsplit) ;

   getint(ph, "minlen:", &minlen) ;
   getint(ph, "maxlen:", &maxlen) ;
   getint(ph, "minlength:", &minlen) ;
   getint(ph, "maxlength:", &maxlen) ;
   getint(ph, "basewindow:", &basewindow) ;
   getint(ph, "mapreg:", &mapreg) ;
   getstring(ph, "cadjustname:", &cadjustname) ;
   getstring(ph, "cadjname:", &cadjustname) ;
   getlongstring(ph, "samxtra:", &samxtra) ;
   getstring(ph, "forcereg:", &forcereg) ; // all reads forced to use this reg.  For prokaryotes
   getint(ph, "killsnps:", &killsnps) ; // kill reference on snps
   getstring(ph, "snpname:", &snpname) ; // needed if killsnps set YES
  

   writeparsx(ph) ;
   closepars(ph) ;
   fflush(stdout) ;

}

int readfa1(char *faname, char **pfasta, int *flen) 
{

 faidx_t *fai ;
 int k, len, t ;
 char *ttfasta, *qqfasta ;
 char ssreg[20] ;
 int ntry = 0, itry ;
 int ismt ;
 
 if (pfasta != NULL) freestring(pfasta) ;

 *flen = 0 ;
 *pfasta = NULL ;

  if (faname  == NULL) { 
   return  0;
  }
  fai = fai_load(faname) ;
  strcpy(ssreg, regname) ;

  ismt = NO ;
  t = strcmp(regname, xregname) ; if ((t==0)  && (mapreg)) strcpy(ssreg, "X") ;
  t = strcmp(regname, yregname) ; if ((t==0) && (mapreg))  strcpy(ssreg, "Y") ;
  t = strcmp(regname, mtregname) ; if ((t==0) && (mapreg))  strcpy(ssreg, "MT") ;
  if (t==0) ismt = YES ; 

  ttfasta = myfai_fetch(fai, ssreg, &len) ;
  *flen = len ;
  if (ttfasta == NULL) { 
   printf ("***warning: no hetfa for %s :: %s\n", faname, ssreg) ;
   return len ;
  }
  if (ismt) { 
// pretend that MT is linear of twice length
   ZALLOC(qqfasta, 2*len+1, char) ; 
   strcpy(qqfasta, ttfasta) ;
   strcpy(qqfasta+len, ttfasta) ;
   *pfasta = strdup(qqfasta) ; 
   len = 2*len ; 
   freestring(&qqfasta) ;
  }
  else {
   *pfasta = strdup(ttfasta) ;
  }

  freestring(&ttfasta) ;
  return len ;

}

char *myfai_fetch(faidx_t *fai, char *reg, int  *plen)
{
  char *treg, *s ;
  treg = strdup(reg) ;

  if (fai==NULL) fatalx("(my_fai_fetch): fai NULL\n") ;

  s = fai_fetch(fai, treg, plen) ;
  if (*plen > 0) {
    free(treg) ;
    return s ;
  }
  free(treg) ;
  return NULL ;
}

void clearfainfo(FATYPE *fapt, int mode)
{

 if (mode==1) {
   freestring(&fapt->faname) ;
   freestring(&fapt->alias) ;
   freestring(&fapt->famask) ;
   fapt -> fai = fapt -> faimask =  NULL ;
 }
 freestring(&fapt -> rstring) ;
 freestring(&fapt -> mstring) ;
 freestring(&fapt -> regname) ;
 fapt -> lopos = fapt -> hipos = -1 ;
 fapt -> len = 0 ;
 fapt -> mlen = 0 ;
 fapt -> rlen = 0 ;

}

void printfapt(FATYPE *fapt)
{
  int k, np = 0, mlen = fapt -> mlen   ;
  char cc ;

  printf("fapt: %s %s\n",  fapt->faname, fapt->alias) ;
  printf("%x %d %d n", fapt -> fai, fapt -> lopos, fapt -> hipos) ;
  printf("len: %d rlen: %d\n", fapt -> len, fapt -> rlen) ;
  printf("mask: %s %d\n", fapt -> famask, fapt -> mlen) ;
  printnl() ;

  fflush(stdout) ;
}
           
int  mksamplist(char **samplist, Indiv **indivmarkers, int numindivs) 
{
  int k, n=0 ;
  Indiv *indx ;

  for (k=0; k<numindivs; ++k) { 
    indx = indivmarkers[k] ;
    if (indx -> ignore) continue ;
    samplist[n] = strdup(indx -> ID) ;
    ++n ;
  }
  return n ;

}

int grabreference(char *reg, int pos, int len, char *refstring)  
{

   char sss[1000] ;
   int t, k, lenx ;
       

   t = 1 ;
   if (regname != NULL) t = strcmpreg(reg, regname) ;
   if (t!=0) {
    loadreg(reg) ; 
    printf("(grabref) loadreg called\n") ;
   }
   lenx = MIN(len, reglen-pos+1) ; 
   if (lenx>999) fatalx("bad length\n") ;
   if (lenx <= 0) { 
    cclear(refstring, '?', len) ; 
    refstring[len] = CNULL ;
    return -1 ;
   }
// last char is regstring + pos -1 + lenx -1 
   if (pos <= 0) return -1 ;
   if ((pos + lenx -2) > (reglen -1)) return -1 ; 
   strncpy(sss, regstring + pos-1, lenx) ;
   for (k=0; k<len-lenx; ++k) { 
    sss[k+lenx] ='?' ;
   }
   sss[len] = CNULL ;
   strcpy(refstring, sss) ;
   return 1 ;
}


void loadreg(char *reg) 
{
  int treglen ;
  int aa[4], t, k, tt ;
  char *faname ;
  char *rreg ; 
  int px[3] ; 

   rreg = strdup(reg) ;

   freestring(&regname) ; 
   freestring(&regstring) ; 
   freestring(&chimpstring) ; 
 
   regname = strdup(rreg) ;
   faname = ref[0] ;

   readfa1(faname, &regstring, &reglen) ;
   printf("+++ reg %s loaded;  length %d\n", rreg, reglen) ;
   if (reglen <= 0) fatalx(" %s fails to load\n", rreg) ;

   if (polymask) {
    faname = polyref[0] ;
    readfa1(faname, &polystring, &polylen) ;
    printf("+++ polymask %s loaded;  length %d\n", rreg, polylen) ;
    if (polylen <= 0) fatalx(" %s fails to load\n", rreg) ;
    reglen = MIN(reglen, polylen) ;
    ivzero(px, 3) ; 
    for (k=0; k< reglen; ++k) {  
     t = 2 ; 
     if (polystring[k] == '1') t = 1 ; 
     if (polystring[k] == '0') t = 0 ; 
     if (t==1)  regstring[k] = '?' ; 
     ++px[t] ; 
//   tt = ranmod(1000) ; 
//   if (tt==0) printf("zzk: %d %d %c %c\n", k, t, regstring[k], polystring[k]) ;
    }
//  printf("zzpx: ") ;  printimat(px, 1, 3) ;
   }

   free(rreg) ;

}

void mkpop (char **rawpops, char **pops, int n) 
{
  char *spt[MAXFF], *sx  ;
  int k, nsplit ;
 
  for (k=0; k<n; ++k) { 
    freestring(&pops[k]) ;
    nsplit = splitupx(rawpops[k], spt, MAXFF, ':') ;
    sx = spt[0] ;
    pops[k] = strdup(sx) ; 
    freeup(spt, nsplit) ;
  }
}

int callsam() 
{
 char sss[MAXSSS] ;
 char fff[MAXSSS] ; 
 char rrr[MAXSSS] ; 
 int t, k, n, nr=0, kk ;
 char *ff ;

 if (regname == NULL) fatalx("regname:  compulsory\n") ;

 sprintf(rrr, "%s", regname) ;
 t = strcmp(regname, xregname) ; if (t==0) strcpy(rrr, "X") ;
 t = strcmp(regname, yregname) ; if (t==0) strcpy(rrr, "Y") ;
 t = strcmp(regname, mtregname) ; if (t==0) strcpy(rrr, "MT") ;
 if (forcereg != NULL) strcpy(rrr, forcereg) ; 
  
 if (bamname == NULL) fatalx("no bamname\n") ;

   if (samxtra == NULL) {
    sprintf(sss, "%s view %s %s ", samtoolspath, bamname, rrr) ;
   }
   else { 
    sprintf(sss, "%s view \"%s\" %s %s ", samtoolspath, samxtra, bamname, rrr) ;
   }
   n = 1 ;

  if (forcereg != NULL) printf("samtools call: %s\n", sss) ;
 fflush(stdout) ;
 fprintf(stderr, "%s\n", sss) ;
 samstream = popen(sss, "r") ;
 printf("samstream opened\n") ; fflush(stdout) ; 
 return n ;
}

int countsam(int firstpos, int lastpos)   
{
  char line[MAXSSS+1] ;
  static int n=0 ;
  int nprint = 0, ret, nraw = 0  ;
  long nn = 0 ; 
  int maxprint = 100 ;
  int ipop, refmatch, k, pp, t, tt, s, rr, x ;
  int pos, ntri = 0, nread = 0   ;
  int a, b ; 
  char ca, cb ;
  char cseq ;
  SNP *cupt ;
  static long nxread = 0, nrgood = 0  ; 
  double ypmd ; 
  static int nbad = 0  ;
  int div, len, zret ;
  int isdel, oldisdel ;

    for (;;) {

     readptx = readlist[1] ; 
     oldreadpt = readlist[2] ; 
     readpt -> parent =  readpt ;     // dangerous bend parent is (void) * ;
     if ((maxread >0) && (nn >= maxread)) {
      printf("maxread limit reached;  rest of data ignored!\n") ; 
      break ;
     }
      else ret = getreadfnorg(readpt, samstream) ;
      ++nxread ; 
      verbose = NO ;
      ++nn ;
      readpt -> mapq = MIN(readpt -> mapq, mqclip) ;

     if (ret == -99) break ;
     if (ret == -98) break ;

//   if (nxread > maxread) enuf("%ld reads processed!\n",  nxread) ;
     if (printrd) printread(readpt) ;


     if (checkcigar(readpt -> cigar) == NO) {  
      ret = readpt -> rckode = -14 ;
      ++badcigar ; 
      if (badcigar <= pubcigar) {
       printf("badcigar: %s\n", readpt -> cigar) ;
       printf("ref: %s\n", readpt -> refseq) ;
       printf("seq: %s\n", readpt -> seq) ;
      }
     }


     oldreadpt = readpt ; 

  
     oldisdel = NO ; 
     if (strlen(readpt -> cigar) > 4) oldisdel = YES ;
    
     isdel = isindel(readpt -> cigar)  ; 
     if (isdel == YES) { 
      ++xindels[1] ; 
     }
     if (isdel == NO) {
      ++xindels[0] ;
     }

     if (ultima && isdel) { 
       zret = fixindels(readpt, readptx) ; 
//     if (isdel != oldisdel) printf("zzqqdel %s %d\n", readpt -> cigar, zret) ;
       if (zret<0) ret = readpt -> rckode = -14 ;
       readptx -> parent = readpt ;  // 
       readpt = readptx ;
       readpt -> isfix = YES ;
       if (ret < 0) { 
        ++nreject ;
        continue ;
       }
     } 


     if (isdel) fixbase(readpt, basewindow) ;
     if (printrd) {
       printf("+++***+++\n") ;
       printread(readptx) ;
       printread(oldreadpt) ;
     } 
   
     ivclip(readpt -> iqual, readpt -> iqual, -999, bqclip, readpt -> readlen) ;

      ++nraw ;

     if (ret == -33) { 
      ++rgmiss ;
     }
     if (readpt -> pos < firstpos) continue ;
// trap ret = -33 
         
     if (ret<0) { 
      ++nreject ;
      continue ;
     }
     pos = readpt -> pos ;
     if (pos > lastpos) break ;

     if (flipstrand) readpt -> strand = 1 - readpt -> strand ;

 //  readpt -> pmdscore = 9 ;  
     ret = checkread(readpt) ;  // key
     if (ret<0) continue ; 
     if (ret==1) { 
      ++nrgood ; 
     } 
     verbose =  NO ;

     t = 0 ; 
     if (hasrggroup) t = readpt -> rgcat = getrgcat(rgcats, rgclist, numrgc, readpt -> rgstring) ;    
     else t = readpt -> rgcat = 0 ; 
     if (t==-9999) ++rgnone ;
     if (t<0) { 
       ++rgmiss ;
       continue ; 
     }
     ++rgcount[t] ;

     bighist = bighistlist[t] ; 
     if (isdel) {
      bighist = bighistlistu[t] ;
      ++nindel ; 
     }
     else ++nindelx ;
     counthist(readpt, bighist) ;

     s = readpt -> strand ; 

     ++nread ; 
    }

   pclose(samstream) ;

   if (nraw==0) printf("*** no data... database problem? ***\n") ;
   fflush(stdout) ;


   return nn ;
}

int checkread(READ *readpt) 
{
// <0 if read rejected
  int t, strand, ipop, ispair ; 
  char c ;
  double ypmd ;

  if ((processdups == NO) &&(readpt -> isdup))  return -7 ;
  if ((ipop = readpt -> ipop) <0) return -9 ;

  if (readpt -> mapq < pthresh -> lomapthresh) return -3 ;
  if (readpt -> mapq > pthresh -> himapthresh) return -3 ;

  ispair = pthresh -> ispair ;

  if ((ispair == YES) && (readpt -> ispair == NO)) return -12 ; 

  if (readpt -> readlen >   maxlen)   return -2 ;
  if (readpt -> readlen <   minlen)   return -2 ;

  t = strlen(readpt -> cigar) ;          
 if (readpt -> isfix == NO) {
  c = readpt -> cigar[t-1] ;
  if (c != 'M') return -1 ;
  if (t>4) return -1 ;
  if (t<2) return -1 ;
 }

  t = indxstring(bad, nbadrg, readpt -> rgstring) ; 
  if (t>=0) return -11 ;


  if (udgtype < 0) return 1 ;
  ypmd = readpt -> pmdscore ;
  if (usepmdscore && (ypmd > -100)) { 
   if (readpt -> pmdscore <= pmdlo) return -4 ;
   if (readpt -> pmdscore >  pmdhi) return -4 ;
  }
// note that score = pmdlo fails => score = pmdlo = 0  = undamaged

  return 1 ;
}

int stripindivs(Indiv **indivmarkers, int numindivs)   
{
  // squeeze out ignore
  // No genotypes
  int n = 0, g, i, k ;
  int x ; 
  Indiv *indx ;

  // n is index of next unused array element 

  for (k=0; k<numindivs; ++k) {  
    if (indivmarkers[k] -> ignore == YES)  continue ;       // don't store
    if (n==k) {                                             // if no ignored found yet, 
     ++n ;                                                  //   next unused is next element
     continue ;                                             //   and no need to copy
    }

    // copy k -> n 
    indx = indivmarkers[n] ;                     // if kth element is not ignored, put it 
    indivmarkers[n] = indivmarkers[k] ;          //   into next unused element
    indx -> idnum = n ;
    ++n ;
  }

  return n ;

}

void clearhist(HIST *hpt) 
{
/**
typedef struct  {
 int mq ; 
 int bq ; 
 int strand ;
 int  pp ;   // - is 5' end  
 long hist[16] ; 
 long hsum ; 
} HIST ;
*/
 hpt -> mq = hpt -> bq = hpt -> strand = -1 ; 
 hpt -> pp = -9999 ;   
 ivzero(hpt -> hist, 16) ; 
 hpt -> hsum = 0 ;
}

int ppp(int rpos, int len)  
{
    int sgn, x, t ;  

    sgn = 1 ; 
    x = rpos + 1  ; 
    t = len - rpos  ; 
    if (t<x) { 
     x = t ; 
     sgn = -1 ; 
    }
    if (x > readedge) return 99 ;  
    return sgn * x ; 
// negative is 3' end of read   

} 

int cbump(long *xbases, int pos, char *rstring, int rlen, int mode)  
{
  long *zbases  ; 
  char cc ;  
  int t ;

  if (pos <= 0) return -2 ;
  if (pos > rlen) return -2 ;
  zbases = xbases + mode*4 ;
  cc = rstring[pos-1] ; 
  t = base2num(cc) ; 
  if (t<0) return -1 ; 
  ++zbases[t] ; 
  return 1 ;
}

void counthist (READ *readpt, HIST **bighist)  
// main count routine
{
     int a, b, k, strand, pp, nn, p, x, t, w, tt ;  
     char ca, cb, cseq, cpoly ; 
     static long nread = 0 ;
     int mq, bq ;
     HIST *hpt ;
     int rlen ; 
     int ffpos, llpos, midpos ;
     long *xbase ;
     int minmq = 20, minbq = 20 ;
     static long nx = 0 ; 
    

     ++nread ;
     rstring = regstring ; 
     rlen = reglen ; 

     if (usechimp) {
      rstring = chimpstring ;
      rlen = chimplen ; 
     }

     if (stopcount) return ; 
     strand = readpt -> strand ;
     mq = readpt -> mapq ;
     mq = iclip(mq, 0, hd1) ;

     xbase = cbases[readpt -> rgcat][strand] ;
     ++nx ;  
     if (nx<-1) {      
      printf("zzcbase: %d %d\n", readpt -> rgcat, strand) ;  
      printlmat(xbase, 3, 4) ;
     }

     k = 0 ; 
     pp = readpt -> pos + k ;   

     cb = cseq = toupper(readpt -> seq[k]) ;        
     ca = rstring[pp-1] ;
     a = base2num(ca) ;
     b = base2num(cb) ; 
     if ((a>=0) && (b>=0)) ++LC[strand][a][b] ;

     k = readpt -> readlen - 1 ;
     ffpos = readpt -> pos - 1 ;  
     pp = readpt -> pos + k ;   
     llpos = pp ;  

     cb = cseq = toupper(readpt -> seq[k]) ;        
     ca = rstring[pp-1] ;
     a = base2num(ca) ;
     b = base2num(cb) ; 
     if ((a>=0) && (b>=0)) ++RC[strand][a][b] ;

     tt = ranmod(100) ; cpoly = '*' ; 
     cpoly = '?' ; 
     for (k=0; k<readpt -> readlen; ++k) { 
      bq = readpt -> iqual[k] ;  
      bq = iclip(bq, 0, hd2) ;
      p = readpt -> pos + k ;   
      if ((p<=0) || (p>rlen)) continue ;
      t = qcpg1(p, rstring, rlen) ; 
      if (killcpg && (t==YES)) { 
       ++cpgmiss ;
       continue ; 
      }
      if (cpgonly && (t==NO)) { 
       ++cpgmiss ; 
       continue ; 
      }
      pp = ppp(k, readpt -> readlen) ; 
      nn = ptoax[pp] ; 
      if (nn < 0) fatalx("bad pp! %d %d\n", pp, nn) ;
      cb = cseq = toupper(readpt -> seq[k]) ;        
      ca = rstring[p-1] ;
      ca = toupper(ca) ; 
      if (polymask && (p<=polylen)) cpoly = polystring[p-1] ; 
      if (cpoly == '1') continue ; 
      a = base2num(ca) ;
      b = base2num(cb) ; 
      if (a<0) continue ;
      if (b<0) continue ;
      mq = MIN(mq, hd1-1) ;
      bq = MIN(bq, hd2-1) ;
      t = htable[mq][bq][strand][nn] ; 
      if (t<0) fatalx("read overlow (htable)  : %d %d %d %d\n", mq, bq, strand, pp) ; 
      hpt = bighist[t] ;  
      w = 4*a + b ;  
      ++hpt -> hist[w] ; 
      ++hpt -> hsum ; 
      if (hpt -> hsum == BIGINT ) stopcount = YES ;  
      maxmq = MAX(maxmq, mq) ;
      maxbq = MAX(maxbq, bq) ;
      if (mq < minmq) continue ; 
      if (bq < minbq) continue ; 
      if (pp==1) cbump(xbase, ffpos, rstring, rlen, 0) ;
      if (pp==-1) cbump(xbase, llpos, rstring, rlen, 1) ;
      if (pp==99) { 
       midpos = p ; 
       cbump(xbase, midpos, rstring, rlen, 2) ;
      }
     }
     if (nx<-1) {      
      printf("zzcbase2: %d %d\n", readpt -> rgcat, strand) ;  
      printlmat(xbase, 3, 4) ;
     }
     
}

void 
clear4D(int ****table, int x1, int x2, int x3, int x4, int val) 
{
   int a1, a2, a3, a4 ;

   for (a1=0; a1<x1; ++a1) { 
    for (a2=0; a2<x2; ++a2) { 
     for (a3=0; a3<x3; ++a3) { 
      for (a4=0; a4<x4; ++a4) { 

     table[a1][a2][a3][a4] = val ;


   }}}}

} 
void csq(long *xa, long *xb) 
{ 

  double zz[8], ychi, ytail ;

  printlmat(xa, 1, 4) ;
  printlmat(xb, 1, 4) ;

  floatitl(zz, xa, 4) ;
  floatitl(zz+4, xb, 4) ;
  if (testconchi(zz, 2, 4) == NO) { 
   printf("*** (csq) warning. No data at read ends\n") ;
   return ;
  }  
  ychi = conchi(zz, 2, 4) ;
  if (ychi < -99) return ; 
  ytail = rtlchsq(1*3, ychi) ; 
  printf("chisq: %9.3f p-value: %12.6f\n", ychi, ytail) ;

}

void cstats(long **cc) 
{
  int npos, t ;
  double ff[4] ; 
  int pp, x ;
  long cca[4], ccb[4], cc99[4]  ;  
  double xc[8] ; 

  npos = 2*readedge + 1 ;

  lvzero(cca, 4) ;
  lvzero(ccb, 4) ;
  lvzero(cc99, 4) ;
  
  for (x=0; x<npos; ++x) {  
   floatitl(ff, cc[x], 4) ;
   if (asum(ff, 4) < 0.1) continue ;
   pp = atop(x, readedge) ;
   printf("cstat %3d :: ", pp) ; 
   printlmatx(cc[x], 1, 4) ; 
   bal1(ff, 4) ;
   printmat(ff, 1, 4) ;
   if (pp==99) lvvp(cc99, cc99, cc[x], 4) ;
   t = abs(pp) ; 
   if ((t>=2) && (t<=4) && (pp>0)) lvvp(cca, cca, cc[x], 4) ;
   if ((t>=2) && (t<=4) && (pp<0)) lvvp(ccb, ccb, cc[x], 4) ;
  }
  printnl() ;
  printf("5' chisq (pos 2-4 + 99)\n") ; 
  csq(cca, cc99) ;
  printnl() ;
  printf("3' chisq (pos 2-4 + 99)\n") ; 
  csq(ccb, cc99) ;
  printnl() ;
}

void writecadj(long **cc, long *cb, int jnum, int snum, FILE *cfff) 
{
  int npos, t ;
  int pp, x ;
  double wscale[4], y, ww[4], wbase[4], ff[4] ; 
  
  if (cfff == NULL) return ;

  fprintf(cfff, "##base distribution: ") ;
  fprintf(cfff, "strand %d stack %d ", snum, jnum) ;
  floatitl(wbase, cb, 4) ; 
  bal1(wbase, 4) ;
  printmatwxfile(wbase, 1, 4, 4, cfff) ; 
  printlmatwfile(cb, 1, 4, 4, cfff) ;
  if (flatbase) {
    fprintf(cfff, "##flatbase used\n") ;
    vclear(wbase, 1.0, 4) ; bal1(wbase, 4) ;
  }
  for (x=0; x<npos; ++x) {  
   floatitl(ff, cc[x], 4) ; 
   vsp(ff, ff, 1, 4) ;
   bal1(ff, 4) ;
   vvd(wscale, wbase, ff, 4) ;  // crude scaling
   vmaxmin(wscale, 4, &y, NULL) ;  
   vst(wscale, wscale, 1.0/y, 4) ;
   pp = atop(x, readedge) ;
   fprintf(cfff, "%3d %3d %3d ", pp, jnum, snum) ; 
   printmatfile(wscale, 1, 4, cfff) ;
  }
}


void countbasepos  ()  
{

  int a1, a2, a3, a4, t, k, j, u, s ; 
  HIST *histpt ; 
  int ds, tt ;
  FILE *fff = stdout ;
  int npos ; 
  long **countb,  *** cbase, **cball ; 
  FILE *cadjfile = NULL ;
  int minmq = 20, minbq = 20 ; 

  if (posbasestats == NO) return ;
  if (cadjustname != NULL) openit (cadjustname, &cadjfile, "w") ; 
  npos = 2*readedge + 1 ;

  countb = initarray_2Dlong(npos, 4, 0) ;
  cbase = initarray_3Dlong(2, numrgcats, 4, 0) ;
  cball = initarray_2Dlong(2, 4, 0) ;

  printf("countb allocated\n") ;

  if (numrgc>1) { 
   for (k=0; k<numrgc; ++k) { 
    fprintf(fff, "##rggroup: ") ; 
    fprintf(fff, "%30s ", rgclist[k]) ; 
    fprintf(fff, " %3d", rgcats[k]) ; 
    t = rg2cat[k] ;
    fprintf(fff, ":  %8d", rgcount[t]) ; 
    fprintf(fff, "\n") ;
   }
   for (k=0; k<numrgc; ++k) { 
    printf( "##rggroup: ") ; 
    printf( "%30s ", rgclist[k]) ; 
    printf( " %3d", rgcats[k]) ; 
    t = rgcats[k] ;
    printf(":  %8d", rgcount[t]) ; 
    printf( "\n") ;
   }
  }
  for (j=0; j<numrgcats; ++j) { 
   for (s=0; s<=1; ++s) {
   bighist = bighistlist[j] ;
   for (a1=0; a1<hd1; ++a1) { 
    for (a2=0; a2<hd2; ++a2) { 
     a3 = s ;
     a4 = 2*readedge ; 
     t = htable[a1][a2][a3][a4] ;  
     if (t<0) continue ; 
     histpt = bighist[t] ; 
     for (k=0; k<16; ++k) { 
       u = k % 4 ;  
       cbase[s][j][u] += histpt -> hist[k] ;
       cball[s][u] += histpt -> hist[k] ;
     } 
   }}}} 
  for (j=0; j<numrgcats; ++j) { 
   if (cadjstackbase) break ;
   for (s=0; s<=1; ++s) {
    copylarr(cball[s], cbase[s][j], 4) ; 
   }}
    



// don't process bighistlistu for now
  for (j=0; j<numrgcats; ++j) { 
   for (s=0; s<=1; ++s) {
   printf("hist %d strand %d\n", j, s) ;
   lclear2D(&countb, npos, 4, 0) ; 
   bighist = bighistlist[j] ;
   for (a1=0; a1<hd1; ++a1) { 
    for (a2=0; a2<hd2; ++a2) { 
     a3 = s ;
      for (a4=0; a4<hd4; ++a4) { 
       t = htable[a1][a2][a3][a4] ;  
       if (t<0) continue ; 
       histpt = bighist[t] ; 
       for (k=0; k<16; ++k) { 
        u = k % 4 ;  
        countb[a4][u] += histpt -> hist[k] ;
       } 
   }}}
   cstats(countb) ;      
   writecadj(countb, cbase[s][j], j, s, cadjfile) ;    
 }} 

 free2Dlong(&countb, npos) ; 
 free3Dlong(&cbase, 2, numrgcats) ; 
 free2Dlong(&cball, 2) ; 

}

void setww(double *a, long *z)  
{
 floatitl(a, z, 12) ;
 bal1(a, 4) ;
 bal1(a+4, 4) ;
 bal1(a+8, 4) ;
}

void dumphist(char *fname, HIST ***bigh, int ultmode) 
{

  FILE *fff ; 
  int a1, a2, a3, a4, t, k, j ; 
  HIST *histpt ; 
  int ds, tt ;
  char fname2[MAXSTR] ;
  static int udgtype = -1 ;
  int numwrite = 0 ; 
  double ww[12] ;

  if (fname == NULL) return ; 

  openit(fname, &fff, "w") ;
  fprintf(fff, "###globals %d %d %d %d", maxmq+1, maxbq+1, hd3, readedge) ; 
  fprintf(fff, " :: %s :: ", sampname) ;
  fprintf(fff, "  %s", bamname) ;
  if (numrgcats >  1) fprintf(fff, " :: %d", numrgcats) ;
  fprintf(fff, "\n") ;

  if (numrgc>1) { 
   for (k=0; k<numrgc; ++k) { 
    fprintf(fff, "##rggroup: ") ; 
    fprintf(fff, "%30s ", rgclist[k]) ; 
    fprintf(fff, " %3d", rgcats[k]) ; 
    t = rgcats[k] ;
    fprintf(fff, ":  %8d", rgcount[t]) ; 
    fprintf(fff, "\n") ;
   }
   for (k=0; k<numrgc; ++k) { 
    printf( "###rggroup: ") ; 
    printf( "%30s ", rgclist[k]) ; 
    printf( " %3d", rgcats[k]) ; 
    t = rgcats[k] ;
    printf( ":  %8d", rgcount[t]) ; 
    printf( "\n") ;
   }
  }

//  printf("zzd: %d %d %s\n", calcudgtype, numrgcats, deaminhistname) ;
  for (j=0; j<numrgcats; ++j) { 
   bighist = bigh[j] ;
   if ((calcudgtype) && (deaminhistname != NULL)) { 
    fprintf(deaminfile, "## deamination probabilities for file %s histogram %d ", fname, j) ; 
    deaminhist(bighist, deaminfile) ;
   }

   if (calcudgtype) {
    if (udgverbose) verbose = YES ;
    if (udgtype <0) calcudg(bighist, &ds, &udgtype) ; 
    fprintf(fff, "## histogram %d ", j) ; 
    fprintf(fff, " : %s\n", udgm(ds, udgtype)) ; 
    verbose = NO ;
   }
   for (k=0; k<=1; ++k) { 
    setww(ww, cbases[j][k]) ; 
    fprintf(fff, "##prevbase") ; printmatwxfile(ww, 1, 4, 4, fff) ;   fprintf(fff, "%d %2d\n", k, j) ;
    fprintf(fff, "##postbase") ; printmatwxfile(ww+4, 1, 4, 4, fff) ; fprintf(fff, "%d %2d\n", k, j) ;
    fprintf(fff, "##midbase ") ; printmatwxfile(ww+8, 1, 4, 4, fff) ; fprintf(fff, "%d %2d\n", k, j) ;
   }

   for (a1=0; a1<=maxmq; ++a1) { 
    for (a2=0; a2<=maxbq; ++a2) { 
     for (a3=0; a3<hd3; ++a3) { 
      for (a4=0; a4<hd4; ++a4) { 
       t = htable[a1][a2][a3][a4] ;  
       if (t<0) continue ; 
       histpt = bighist[t] ; 
       if (histpt -> hsum == 0) continue ; 
       fprintf(fff, "%2d %2d %d %3d ", histpt -> mq, histpt -> bq, histpt -> strand, histpt -> pp) ; 
       for (k=0; k<16; ++k) { 
        fprintf(fff, " %6ld", histpt -> hist[k]) ;
        t = k % 4 ;  if ((t == 3) && (k<15)) fprintf(fff, "  ") ; 
       } 
       if (numrgcats >  1) fprintf(fff, " :: %d", j) ;
       fprintf(fff, "\n") ;
       ++numwrite ;
   }}}}
 }
  fclose(fff) ;
  printf("## histfile %s %d data lines written\n", fname, numwrite) ;


  if (ultmode == 1) { 
   mkfname2(fname2, fname) ; 
   dumphist(fname2, bighistlistu, 0) ;
  }

}

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

 HIST *histpt  ; 
 int mapq, baseq,strand, pos, ppos ;   
 int a1, a2, a3, a4 ; 
  int j, jmax, t, dlen ;
  long *bbb, *ttt  ; 
  long **bdata, ***bbdata ; 
  int pp, x ;
  int *x2p ;

  double aa[4], yct, yga, ychi, yrat ; 
  double yca, ycb ; 
  int ds = NO, ph = YES, tt = 1  ; 
  double y1, y2, *ww ; 
  double **deamin, dabase[2] ;
  long bb[8] ; 

 
 dlen = 2*readedge + 1 ; 

 ZALLOC(x2p, dlen+5, int) ;

 ZALLOC(bbdata, 2, long **) ; 
 bbdata[0] = initarray_2Dlong(dlen+5, 4, 0) ;
 bbdata[1] = initarray_2Dlong(dlen+5, 4, 0) ;

 deamin = initarray_2Ddouble(dlen+1, 2, 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) { 
       t = htable[a1][a2][a3][a4] ;  
       if (t<0) continue ; 
       histpt = bighist[t] ; 
       if (histpt -> hsum == 0) continue ; 
       crackx(bb, histpt) ;
       strand = bb[2] ;  
       pp = bb[3] ;  
       x = ptox(pp) ; 
       x2p[x] = pp ; 
       ttt = bbdata[strand][x] ; 
       lvvp(ttt, ttt, bb+4, 4) ;
   }}}}

 verbose = NO ;
 fprintf(fff, " (neg pos => 3' end) \n") ;
// fprintf(fff, "%8s %3s ", "strand", "pos") ;  
// fprintf(fff, "Naive          Adjusted\n") ;
 fprintf(fff, "strand pos         Naive                Adjusted\n") ;
 fprintf(fff, "               C -> T    G -> A       C -> T   G -> A") ;
 fprintf(fff, "\n") ;
// fprintf(fff, "zz %d\n", dlen) ;
 for (strand=0; strand <=1; ++strand) {
  for (x=1; x<=dlen; ++x) {  
   pp = x2p[x] ; 
   if (verbose) { 
    printf(" %3d %3d : ", strand, pp) ;
    printlmat(bbdata[strand][x], 1, 4) ;
   } 
   floatitl(aa, bbdata[strand][x], 4)  ; 
   vsp(aa, aa, 1, 4) ; 
   yct = aa[1]/asum(aa, 2) ; 
   yga = aa[3]/asum(aa+2, 2) ; 
   deamin[x][0] = yct ; 
   deamin[x][1] = yga ; 
   if (pp == 99) copyarr(deamin[x], dabase, 2) ;
  }
  for (x=1; x<=dlen; ++x) { 
   pp = x2p[x] ; 
//   if (pp==99) continue ;
   ww = deamin[x] ; 
   fprintf(fff, "%8d %3d ", strand, pp) ;
   fprintf(fff, " %12.6f", ww[0]) ; 
   fprintf(fff, " %12.6f", ww[1]) ; 
   fprintf(fff, "  ") ; 
   vvm(ww, ww, dabase, 2) ; 
   vclip(ww, ww, 0, 999, 2) ;
   fprintf(fff, " %12.6f", ww[0]) ; 
   fprintf(fff, " %12.6f", ww[1]) ; 
   fprintf(fff, "\n") ;
  }
 }

 verbose = NO ;

 free2Dlong(&bbdata[0], dlen+5) ;
 free2Dlong(&bbdata[1], dlen+5) ;
 free(bbdata) ;
 free2D(&deamin, dlen+1) ;
 free(x2p) ;


 return  ;



}

/**
void calcudg(HIST **bighist, int *stype, int *utype) 
// calculate type for a histogram 
{

#define MAXCL 20000 
 HIST *histpt  ; 
 int mapq, baseq,strand, pos, ppos ;   
 long cc[4], bb[8], **clist, *csum ;  
 int maxcl, numcl ; 
 int a1, a2, a3, a4 ; 
  int j, jmax, t, dlen ;
  long *bbb, *ttt  ; 
  long **bdata, ***bbdata ; 
  int pp, x ;
  int *x2p ;

  double aa[4], yct, yga, ychi, yrat ; 
  double yca, ycb ; 
  char *strandtype, *udgmode ;
  int ds = NO, ph = YES, tt = 1  ; 
  double y1, y2 ; 


 maxcl = hd1*hd2*hd3*hd4 + 1000 ; 
 printf("calcudg verbose: %d\n") ;
 
 dlen = 2*readedge + 1 ; 

 ZALLOC(bbdata, 2, long **) ; 
 bbdata[0] = initarray_2Dlong(dlen+5, 4, 0) ;
 bbdata[1] = initarray_2Dlong(dlen+5, 4, 0) ;

  numcl = 0 ; 

   clist = initarray_2Dlong(maxcl, 8, 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) { 
       t = htable[a1][a2][a3][a4] ;  
       if (t<0) continue ; 
       histpt = bighist[t] ; 
       if (histpt -> hsum == 0) continue ; 
       crackx(bb, histpt) ;
       copylarr(bb, clist[numcl], 8) ;
        ++numcl ;
   }}}}

 ZALLOC(csum, numcl, long) ; 
 for (j=0; j<numcl; ++j) { 
  csum[j] = longsum(clist[j]+4, 4) ; 
 }
 lvlmaxmin(csum, numcl, &jmax, NULL) ;
 bbb = clist[jmax] ;  
 mapq = bbb[0] ; 
 baseq = bbb[1] ;  
 strand = bbb[2] ;  
 pp = bbb[3] ;  

 ZALLOC(x2p, dlen+5, int) ; 
 for (j=0; j<numcl; ++j) { 
   ttt = clist[j] ; 
   if (ttt[0] != mapq) continue ;
   if (ttt[1] != baseq) continue ;
   strand = ttt[2] ; 
   if (strand<0) fatalx("bad strand\n")  ; 
   if (strand>1) fatalx("bad strand\n")  ; 
   pp = ttt[3]; 
   x = ptox(pp) ;
   x2p[x] = pp ;
   copylarr(ttt+4, bbdata[strand][x], 4) ; 
 }
 if (verbose) { 
  printf("%10s", "##") ;
  printf("%10s ", "CC") ; 
  printf("%10s ", "CT") ; 
  printf("%10s ", "GG") ; 
  printf("%10s ", "GA") ; 
  printnl() ;
 } 
 
 strandtype = strdup("single") ; 
 udgmode = strdup("half") ; 
 tt = 1 ; // udg 
 for (strand=0; strand <=1; ++strand) {
  for (x=1; x<=dlen; ++x) {  
   pp = x2p[x] ; 
   if (verbose) { 
    printf(" %3d %3d : ", strand, pp) ;
    printlmatx(bbdata[strand][x], 1, 4) ;
   } 

   floatitl(aa, bbdata[strand][x], 4)  ; 
   vsp(aa, aa, 1, 4) ; 
   yct = aa[1]/asum(aa, 2) ; 
   yga = aa[3]/asum(aa+2, 2) ; 
   if ((strand == 0) && (pp == 1)) yca = yct ; 
   if ((strand == 0) && (pp == 5)) ycb = yct ; 
   if (verbose) {
    printf(" %9.3f", 1000*yct) ;
    printf(" %9.3f", 1000*yga) ;
   } 
   yrat = log(yct/yga)  ;
   ychi = conchi(aa, 2, 2) ;     
   if (verbose) {
    printf(" %9.3f", yrat) ; 
    printf(" %9.3f", ychi) ; 
    printnl() ;
   } 
   if ((pp==99) && (fabs(yrat)<0.5)) {
     strandtype = strdup("double") ;
     ds = YES ; 
   } 
   if ((strand == 0) && (pp==5) && (ycb>0.1*yca)) {
    udgmode = strdup("minus") ;
    ph = NO ;
    tt = 0 ; 
   }
  }
 } 
 if (tt>0) { 
  strand = 0 ; x = 2 ; floatitl(aa, bbdata[strand][x], 4)  ;  vsp(aa, aa, 1, 4) ; y2 = aa[1]/asum(aa, 2) ;
  strand = 0 ; x = 1 ; floatitl(aa, bbdata[strand][x], 4)  ;  vsp(aa, aa, 1, 4) ; y1 = aa[1]/asum(aa, 2) ;
  yrat = y1/y2 ; 
  if (yrat < 4.0) tt = 2 ; 
 }
 printf("## strand udg: %s\n", udgm(ds, tt)) ; 

 free2Dlong(&clist, maxcl) ;
 free(csum) ;
 free2Dlong(&bbdata[0], dlen+5) ;
 free2Dlong(&bbdata[1], dlen+5) ;
 free(bbdata) ;
 free(x2p) ;

 *utype = tt ; 
 *stype = ds ; 

 return  ;
}
*/

char *udgm(int ds, int tt) 
{
 char ss[20], *s1, *s2 ; 
 if (ds == YES) s1 = strdup("double") ;
 if (ds == NO) s1 = strdup("single") ;
 if (tt == 0) s2 = strdup("minus") ;
 if (tt == 1) s2 = strdup("half") ;
 if (tt == 2) s2 = strdup("plus") ;

 sprintf(ss, "%s:%s", s1, s2) ; 
 return strdup(ss) ; 


}

int ptox(int p) 
{ 
 int t ;  
 int dlen ;

 dlen = 2*readedge + 1 ;
  
 if (p == 99) return dlen ;
 if (p<0) { 
  t = -p ; 
  return dlen - t  ;
 }
 return p  ;  
}
void crackx(long *bb, HIST *hpt )  
{ 

 long  *bx ;
 int *jx ; 
 int k ;  

 jx = hpt -> hist ;  
 bb[0] = hpt -> mq ; 
 bb[1] = hpt -> bq ; 
 bb[2] = hpt -> strand ; 
 bb[3] = hpt -> pp ; 
 bx = bb + 4 ; 
 bx[0] = jx[5] ;  // CC 
 bx[1] = jx[7] ;  // CT 
 bx[2] = jx[10];  // GG 
 bx[3] = jx[8] ;  // GA 

}

void calcudg(HIST **bighist, int *stype, int *utype) 
// calculate type for a histogram 
{

#define MAXCL 20000 
 HIST *histpt  ; 
 int mapq, baseq,strand, pos, ppos ;   
 long cc[4], bb[8], **clist, *csum ;  
 int maxcl, numcl ; 
 int a1, a2, a3, a4 ; 
  int j, jmax, t, dlen ;
  long *bbb, *ttt  ; 
  long **bdata, ***bbdata ; 
  int pp, x ;
  int *x2p ;

  double aa[4], yct, yga, ychi, yrat ; 
  double yca, ycb ; 
  char *strandtype, *udgmode ;
  int ds = NO, ph = YES, tt = 1  ; 
  double y1, y2 ; 


 maxcl = hd1*hd2*hd3*hd4 + 1000 ; 
 
 dlen = 2*readedge + 1 ; 

 ZALLOC(bbdata, 2, long **) ; 
 bbdata[0] = initarray_2Dlong(dlen+5, 4, 0) ;
 bbdata[1] = initarray_2Dlong(dlen+5, 4, 0) ;

  numcl = 0 ; 

   clist = initarray_2Dlong(maxcl, 8, 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) { 
       t = htable[a1][a2][a3][a4] ;  
       if (t<0) continue ; 
       histpt = bighist[t] ; 
       if (histpt -> hsum == 0) continue ; 
       crackx(bb, histpt) ;
       copylarr(bb, clist[numcl], 8) ;
        ++numcl ;
   }}}}

 ZALLOC(csum, numcl, long) ; 
 for (j=0; j<numcl; ++j) { 
  csum[j] = longsum(clist[j]+4, 4) ; 
 }
 lvlmaxmin(csum, numcl, &jmax, NULL) ;
 bbb = clist[jmax] ;  
 mapq = bbb[0] ; 
 baseq = bbb[1] ;  
 strand = bbb[2] ;  
 pp = bbb[3] ;  

 if (verbose) { 
  printf("calcudg: largest histogram line ") ; printlmat(bbb, 1, 4) ;
 }

 ZALLOC(x2p, dlen+5, int) ; 
 for (j=0; j<numcl; ++j) { 
   ttt = clist[j] ; 
   if (ttt[0] != mapq) continue ;
   if (ttt[1] != baseq) continue ;
   strand = ttt[2] ; 
   if (strand<0) fatalx("bad strand\n")  ; 
   if (strand>1) fatalx("bad strand\n")  ; 
   pp = ttt[3]; 
   x = ptox(pp) ;
   x2p[x] = pp ;
   copylarr(ttt+4, bbdata[strand][x], 4) ; 
 }
 strandtype = strdup("single") ; 
 udgmode = strdup("half") ; 
 tt = 1 ; // udg  half
 for (strand=0; strand <=1; ++strand) {
  for (x=1; x<=dlen; ++x) {  
   pp = x2p[x] ; 
   if (verbose) { 
    printf(" %3d %3d : ", strand, pp) ;
    printlmatx(bbdata[strand][x], 1, 4) ;
   } 

   floatitl(aa, bbdata[strand][x], 4)  ; 
   vsp(aa, aa, 1, 4) ; 
   yct = aa[1]/asum(aa, 2) ; 
   yga = aa[3]/asum(aa+2, 2) ; 
   if ((strand == 0) && (pp == 1)) yca = yct ; 
   if ((strand == 0) && (pp == 5)) ycb = yct ; 
   if (verbose) {
    printf(" %9.3f", 1000*yct) ;
    printf(" %9.3f", 1000*yga) ;
   } 
   yrat = log(yct/yga)  ;
   ychi = conchi(aa, 2, 2) ;     
   if (verbose) {
    printf(" yca %9.3f", 1000*yca) ; 
    printf(" ycb %9.3f", 1000*ycb) ; 
    printf(" yrat %9.3f", yrat) ; 
    printf(" ychi %9.3f", ychi) ; 
    printnl() ;
   } 
   if ((pp==99) && (fabs(yrat)<0.5)) {
     strandtype = strdup("double") ;
     ds = YES ; 
   } 
   if ((strand == 0) && (pp==5) && (ycb>0.1*yca)) {
    udgmode = strdup("minus") ;
    ph = NO ;
    tt = 0 ; 
   }
  }
 } 
 if (tt>0) { 
  strand = 0 ; x = 2 ; floatitl(aa, bbdata[strand][x], 4)  ;  vsp(aa, aa, 1, 4) ; y2 = aa[1]/asum(aa, 2) ;
  strand = 0 ; x = 1 ; floatitl(aa, bbdata[strand][x], 4)  ;  vsp(aa, aa, 1, 4) ; y1 = aa[1]/asum(aa, 2) ;
  yrat = y1/y2 ; 
  if (yrat < 4.0) tt = 2 ; 
 }
 printf("## strand udg: %s\n", udgm(ds, tt)) ; 

 free2Dlong(&clist, maxcl) ;
 free(csum) ;
 free2Dlong(&bbdata[0], dlen+5) ;
 free2Dlong(&bbdata[1], dlen+5) ;
 free(bbdata) ;
 free(x2p) ;

 *utype = tt ; 
 *stype = ds ; 

 return  ;
}