pbwtIO.cpp

/*  File: pbwtIO.c
 *  Author: Richard Durbin (rd@sanger.ac.uk)
 *  Copyright (C) Genome Research Limited, 2013-
 *-------------------------------------------------------------------
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *   http://www.apache.org/licenses/LICENSE-2.0
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *-------------------------------------------------------------------
 * Description: read/write functions for pbwt package
 * Exported functions:
 * HISTORY:
 * Last edited: Mar 11 16:38 2016 (rd)
 * readPhase updated for chromopainter and chromopainter v2 formats
 * Created: Thu Apr  4 11:42:08 2013 (rd)
 *-------------------------------------------------------------------
 */

#include "pbwt.h"
#include <ctype.h>

int nCheckPoint = 0 ;	/* if set non-zero write pbwt and sites files every n sites when parsing external files */

static BOOL isWriteImputeRef = FALSE ;	/* modifies WriteSites() and WriteHaplotypes() for pbwtWriteImputeRef */

/* basic function to store packed PBWT */

void pbwtWrite (PBWT *p, FILE *fp) /* just writes compressed pbwt in yz */
{
  if (!p || !p->yz) die ("pbwtWrite called without a valid pbwt") ;
  if (!p->aFstart || !p->aFend) die ("pbwtWrite called without start and end indexes") ;
  /* version 2 added start and end indexes */
  if (fwrite ("PBW3", 1, 4, fp) != 4) /* version 3 with 8 byte pbwt size */
    die ("error writing PBWT in pbwtWrite") ;
  if (fwrite (&p->M, sizeof(int), 1, fp) != 1)
    die ("error writing M in pbwtWrite") ;
  if (fwrite (&p->N, sizeof(int), 1, fp) != 1)
    die ("error writing N in pbwtWrite") ;
  if (fwrite (p->aFstart, sizeof(int), p->M, fp) != p->M)
    die ("error writing aFstart in pbwtWrite") ;
  if (fwrite (p->aFend, sizeof(int), p->M, fp) != p->M)
    die ("error writing aFend in pbwtWrite") ;
  long n = arrayMax(p->yz) ;
  if (fwrite (&n, sizeof(long), 1, fp) != 1)
    die ("error writing n in pbwtWrite") ;
  if (fwrite ("    ", 1, 4, fp) != 4)
    die ("error writing padding space in pbwtWrite") ;
  if (fwrite (arrp(p->yz, 0, uchar), sizeof(uchar), arrayMax(p->yz), fp) != arrayMax(p->yz))
    die ("error writing data in pbwtWrite") ;

  fprintf (logFile, "written %ld chars pbwt: M, N are %d, %d\n", arrayMax(p->yz), p->M, p->N) ;
}

void pbwtWriteSites (PBWT *p, FILE *fp)
{
  if (!p || !p->sites) die ("pbwtWriteSites called without sites") ;

  int i ;
  for (i = 0 ; i < p->N ; ++i)
    { Site *s = arrp(p->sites, i, Site) ;
      if (isWriteImputeRef)
	fprintf (fp, "site%d\t%d", i+1, s->x) ;
      else
	fprintf (fp, "%s\t%d", p->chrom ? p->chrom : ".", s->x) ;
      fprintf (fp, "\t%s", dictName (variationDict, s->varD)) ;
      fputc ('\n', fp) ;
    }
  if (ferror (fp)) die ("error writing sites file") ;

  fprintf (logFile, "written %d sites from %d to %d\n", p->N, 
	   arrp(p->sites, 0, Site)->x, arrp(p->sites, p->N-1, Site)->x) ;
}

void pbwtWriteSamples (PBWT *p, FILE *fp)
{
  if (!p || !p->samples) die ("pbwtWriteSamples called without samples") ;

  int i ;
  for (i = 0 ; i < p->M ; i += 2) /* assume diploid for now */
    { Sample *s = sample (p, i) ;
      fprintf (fp, "%s", sampleName(s)) ;
      if (s->popD) fprintf (fp, "\tPOP:%s", popName(s)) ;
      if (s->mother) fprintf (fp, "\tMOTHER:%s", sampleName(sample (p, s->mother))) ;
      if (s->father) fprintf (fp, "\tFATHER:%s", sampleName(sample (p, s->father))) ;
      fputc ('\n', fp) ;
    }     
  if (ferror (fp)) die ("error writing samples file") ;

  fprintf (logFile, "written %d samples\n", p->M/2) ;
}

void writeDataOffset (FILE *fp, char *name, Array data, Array offset, int N)
{
  if (!offset || !data) die ("write %s called without data", name) ;
  int dummy = -1 ;	/* ugly hack to mark that we now use longs not ints */
  if (fwrite (&dummy, sizeof(int), 1, fp) != 1)
    die ("error writing marker in write %s", name) ;
  long n = arrayMax(data) ;
  if (fwrite (&n, sizeof(long), 1, fp) != 1)
    die ("error writing n in write %s", name) ;
  if (fwrite (arrp(data, 0, uchar), sizeof(uchar), n, fp) != n)
    die ("error writing data in write %s", name) ;
  if (fwrite (arrp(offset, 0, long), sizeof(long), N, fp) != N)
    die ("error writing offsets in write %s", name) ;

  fprintf (logFile, "written %ld chars compressed %s data\n", n, name) ;
}

void pbwtWriteMissing (PBWT *p, FILE *fp)
{ writeDataOffset (fp, "missing", p->zMissing, p->missingOffset, p->N) ; }

void pbwtWriteDosage (PBWT *p, FILE *fp)
{ writeDataOffset (fp, "dosage", p->zDosage, p->dosageOffset, p->N) ; }

void pbwtWriteReverse (PBWT *p, FILE *fp)
{
  if (!p || !p->zz) die ("pbwtWriteReverse called without reverse pbwt") ;

  Array tz = p->yz ; p->yz = p->zz ;
  int* tstart = p->aFstart ; p->aFstart = p->aRstart ;
  int* tend = p->aFend ; p->aFend = p->aRend ;

  fprintf (logFile, "reverse: ") ; pbwtWrite (p, fp) ;
  
  p->yz = tz ; p->aFstart = tstart ; p->aFend = tend ;
}


void pbwtWriteAll (PBWT *p, char *root)
{
  FILE *fp ;
#define FOPEN_W(tag)  if (!(fp = fopenTag (root, tag, "w"))) die ("failed to open root.%s", tag)
  FOPEN_W("pbwt") ; pbwtWrite (p, fp) ; fclose (fp) ;
  if (p->sites) { FOPEN_W("sites") ; pbwtWriteSites (p, fp) ; fclose (fp) ; }
  if (p->samples) { FOPEN_W("samples") ; pbwtWriteSamples (p, fp) ; fclose (fp) ; }
  if (p->missingOffset) { FOPEN_W("missing") ; pbwtWriteMissing (p, fp) ; fclose (fp) ; }
  if (p->dosageOffset) { FOPEN_W("dosage") ; pbwtWriteDosage (p, fp) ; fclose (fp) ; }
  if (p->zz) { FOPEN_W("reverse") ; pbwtWriteReverse (p, fp) ; fclose (fp) ; }
}

void pbwtWritePhase (PBWT *p, char *filename)
{
  FILE *fp ;
  if (!(fp = fopen (filename, "w"))) die ("failed to open %s",filename);
  fprintf(fp,"%i\n",p->M);
  fprintf(fp,"%i\nP",p->N);
  int i ; for (i = 0 ; i < p->N ; i++) fprintf(fp," %i",arrayp(p->sites,i,Site)->x);
  fprintf(fp,"\n");

  pbwtWriteTransposedHaplotypes(p,fp);
}

void pbwtCheckPoint (PbwtCursor *u, PBWT *p)
{
  static BOOL isA = TRUE ;
  char fileNameRoot[20] ;

  pbwtCursorToAFend (u, p) ;
  sprintf (fileNameRoot, "check_%c", isA ? 'A' : 'B') ;
  pbwtWriteAll (p, fileNameRoot) ;

  isA = !isA ;
}

/*******************************/

PBWT *pbwtRead (FILE *fp) 
{
  int m, n ;
  long nz ;
  PBWT *p ;
  static char tag[5] = "test" ;
  char pad[4] ;
  int version ;

  if (fread (tag, 1, 4, fp) != 4) die ("failed to read 4 char tag - is file readable?") ;
  if (!strcmp (tag, "PBW3")) version = 3 ; /* current version */
  else if (!strcmp (tag, "PBW2")) version = 2 ; /* with 4 byte count */
  else if (!strcmp (tag, "PBWT")) version = 1 ; /* without start, end indexes */
  else if (!strcmp (tag, "GBWT")) version = 0 ; /* earliest version */
  else die ("failed to recognise file type %s in pbwtRead - was it written by pbwt?", tag) ;

  if (fread (&m, sizeof(int), 1, fp) != 1) die ("error reading m in pbwtRead") ;
  if (fread (&n, sizeof(int), 1, fp) != 1) die ("error reading n in pbwtRead") ;
  p = pbwtCreate (m, n) ;
  if (version > 1)		/* read aFstart and aFend */
    { p->aFstart = myalloc (m, int) ;
      if (fread (p->aFstart, sizeof(int), m, fp) != m) die ("error reading aFstart in pbwtRead") ;
      p->aFend = myalloc (m, int) ;
      if (fread (p->aFend, sizeof(int), m, fp) != m) die ("error reading aFend in pbwtRead") ;
    }
  else				/* set aFstart to 0..M-1, leave aFend empty */
    { p->aFstart = myalloc (m, int) ;
      int i ; for (i = 0 ; i < m ; ++i) p->aFstart[i] = i ;
    }

  if (version <= 2)
    { if (fread (&n, sizeof(int), 1, fp) != 1) die ("error reading pbwt file") ;
      nz = n ;
    }
  else
    if (fread (&nz, sizeof(long), 1, fp) != 1 ||
	fread (pad, 1, 4, fp) != 4) die ("error reading pbwt file") ;

  p->yz = arrayCreate (nz, uchar) ;
  array(p->yz, nz-1, uchar) = 0 ; /* sets arrayMax */
  if (fread (arrp(p->yz, 0, uchar), sizeof(uchar), nz, fp) != nz)
    die ("error reading data in pbwt file") ;

  fprintf (logFile, "read pbwt %s file with %ld bytes: M, N are %d, %d\n", tag, nz, p->M, p->N) ;
  return p ;
}

static BOOL readMatchChrom (char **pChrom, FILE *fp)
{
  char *newChrom = fgetword (fp) ;

  if (strcmp (newChrom, "."))
    { if (!*pChrom) 
	*pChrom = strdup (newChrom) ;
      else if (strcmp (newChrom, *pChrom)) 
	return FALSE ;
    }
  return TRUE ;
}

Array pbwtReadSitesFile (FILE *fp, char **chrom)
{
  char c ;
  Site *s ;
  int line = 1 ;
  Array varTextArray = arrayCreate (256, char) ;
  Array sites = arrayCreate (4096, Site) ;

  while (!feof(fp))
    if (readMatchChrom (chrom, fp))	/* if p->chrom then match, else set if not "." */
      { if (feof(fp)) break ;
	s = arrayp(sites, arrayMax(sites), Site) ;
	s->x = 0 ; while (isdigit(c = fgetc(fp))) s->x = s->x*10 + (c-'0') ;
	if (!feof(fp) && c != '\n')
	  { if (!isspace (c)) die ("bad position line %d in sites file", line) ;
	    while (isspace(c = fgetc(fp)) && c != '\n') ;
	    if (c == '\n') die ("bad end of line at line %d in sites file", line) ;
	    int i = 0 ; array(varTextArray, i++, char) = c ;
	    while ((c = fgetc(fp)) && c != '\n') 
	      array(varTextArray, i++, char) = c ;
	    array(varTextArray, i, char) = 0 ;
	    dictAdd (variationDict, arrayp(varTextArray,0,char), &s->varD) ;
	    while (c != '\n' && !feof(fp)) c = fgetc(fp) ;
	  }
	++line ;
      }
    else if (!feof(fp))
      die ("failed to match chromosome in sites file: line %d", line) ;

  if (ferror (fp)) die ("error reading sites file") ;
  
  fprintf (logFile, "read %ld sites on chromosome %s from file\n", arrayMax(sites), *chrom) ;

  arrayDestroy (varTextArray) ;
  return sites ;
}

void pbwtReadSites (PBWT *p, FILE *fp)
{
  if (!p) die ("pbwtReadSites called without a valid pbwt") ;

  p->sites = pbwtReadSitesFile (fp, &p->chrom) ;
  if (arrayMax(p->sites) != p->N)
    die ("sites file contains %ld sites not %d as in pbwt", arrayMax(p->sites), p->N) ;
}

Array readRefFreq (PBWT *p, FILE *fp)
{
  Array a = arrayCreate (p->N, Site) ;
  char chrom[256], var[256] ;
  int pos ;
  double freq ;
  while (!feof(fp))
    { if (fscanf (fp, "%s\t%d\t%lf\t", chrom, &pos, &freq) != 3 && !feof (fp))
	die ("can't read ref freq file properly") ;
      if (strcmp (chrom, p->chrom))
	die ("chromosome mismatch in readRefFreq '%s' is not '%s'", chrom, p->chrom) ;
      char *cp = var ; while ((*cp = getc(fp)) != '\n' && !feof(fp)) ++cp ; *cp = 0 ;
      Site *s = arrayp(a, arrayMax(a), Site) ;
      s->x = pos ; 
      s->refFreq = freq ;
      dictAdd (variationDict, var, &s->varD) ;
    }
  return a ;
}

void pbwtReadRefFreq (PBWT *p, FILE *fp)
{ if (!p || !p->sites) die ("pbwtReadRefFreq called without current site information") ;
  Array a = readRefFreq (p, fp) ;
  
  int i = 0, j = 0 ;
  Site *ps = arrp(p->sites,i,Site), *as = arrp(a,j,Site) ;
  while (i < p->N)
    { while (j < arrayMax(a) && 
	     (as->x < ps->x || (as->x == ps->x && as->varD < ps->varD))) { ++j ; ++as ; }
      if (ps->x == as->x && ps->varD == as->varD) ps->refFreq = as->refFreq ;
      ++i ; ++ps ;
    }
}

Array pbwtReadSamplesFile (FILE *fp) /* for now assume all samples diploid */
/* should add code to this to read father and mother and population
   propose to use IMPUTE2 format for this */
{
  char *name, c ;
  int line = 0 ;
  Array nameArray = arrayCreate (1024, char) ;
  Array samples = arrayCreate (1024, int) ;

  while (!feof(fp))
    { int n = 0 ;
      while ((c = getc(fp)) && !isspace(c) && !feof (fp)) array(nameArray, n++, char) = c ;
      if (feof(fp)) break ;
      if (!n) die ("no name line %ld in samples file", arrayMax(samples)+1) ;
      array(nameArray, n++, char) = 0 ;
      /* next bit of code deals with header lines for IMPUTE2 samples file, so can read that */
      if (!strcmp (arrp(nameArray,0,char), "ID_1") && !arrayMax(samples))
	{ while ((c = getc(fp)) && c != '\n' && !feof(fp)) ; /* remove header line */
	  while ((c = getc(fp)) && c != '\n' && !feof(fp)) ; /* and next line of zeroes?? */
	  continue ;
	}
      array(samples,arrayMax(samples),int) = sampleAdd (arrp(nameArray,0,char), 0, 0, 0) ;
      /* now remove the rest of the line, for now */
      while (c != '\n' && !feof(fp)) c = getc(fp) ;
    }
  arrayDestroy (nameArray) ;

  fprintf (logFile, "read %ld sample names\n", arrayMax(samples)) ;

  return samples ;
}

void pbwtReadSamples (PBWT *p, FILE *fp)
{
  if (!p) die ("pbwtReadSamples called without a valid pbwt") ;
  Array samples = pbwtReadSamplesFile (fp) ;
  if (arrayMax(samples) != p->M/2) 
    die ("wrong number of diploid samples: %d needed", p->M/2) ;
  p->samples = arrayReCreate(p->samples, p->M, int) ;
  int i ; 
  for (i = 0 ; i < arrayMax(samples) ; ++i)
    { array(p->samples, 2*i, int) = arr(samples, i, int) ;
      array(p->samples, 2*i+1, int) = arr(samples, i, int) ;
    }
  arrayDestroy (samples) ;
}

static void readDataOffset (FILE *fp, char *name, Array *data, Array *offset, int N)
{
  long n ;			/* size of data file */
  int dummy ; 
  if (fread (&dummy, sizeof(int), 1, fp) != 1) 
    die ("read error in read %s", name) ;
  if (dummy != -1) n = dummy ;	/* old version with ints not longs */
  else if (fread (&n, sizeof(long), 1, fp) != 1) 
    die ("read error in read %s", name) ;

  *data = arrayReCreate (*data, n, uchar) ;
  if (fread (arrp(*data, 0, uchar), sizeof(uchar), n, fp) != n)
    die ("error reading z%s in pbwtRead%s", name, name) ;
  arrayMax(*data) = n ;
  fprintf (logFile, "read %ld chars compressed %s data\n", n, name) ;

  *offset = arrayReCreate (*offset, N, long) ;
  if (dummy != -1)		/* old version with ints not longs */
    { /* abuse *offset to hold ints, then update in place */
      if (fread (arrp(*offset, 0, int), sizeof(int), N, fp) != N) 
	die ("error reading %s in pbwtRead%s", name, name) ;
      for (n = N ; n-- ;) 
	arr(*offset, n, long) = arr(*offset, n, int) ; /* !! */
    }
  else
    if (fread (arrp(*offset, 0, long), sizeof(long), N, fp) != N)
      die ("error reading %s in pbwtRead%s", name, name) ;
  arrayMax(*offset) = N ;
}

void pbwtReadMissing (PBWT *p, FILE *fp)
{ readDataOffset (fp, "missing", &p->zMissing, &p->missingOffset, p->N) ; }

void pbwtReadDosage (PBWT *p, FILE *fp)
{ readDataOffset (fp, "dosage", &p->zDosage, &p->dosageOffset, p->N) ; }

void pbwtReadReverse (PBWT *p, FILE *fp)
{
  if (!p) die ("pbwtReadReverse called without a valid pbwt") ;

  PBWT *q = pbwtRead (fp) ;
  if (q->M != p->M || q->N != p->N)
    die ("M %d or N %d in reverse don't match %, %d in forward", q->M, q->N, p->M, p->N) ;
  p->zz = q->yz ; q->yz = 0 ;
  p->aRstart = q->aFstart ; q->aFstart = 0 ;
  p->aRend = q->aFend ; q->aFend = 0 ;
  pbwtDestroy (q) ;
 }

PBWT *pbwtReadAll (char *root)
{
  PBWT *p ;

  FILE *fp ;
  if ((fp = fopenTag (root, "pbwt", "r"))) { p = pbwtRead (fp) ; fclose (fp) ; } 
  else die ("failed to open %s.pbwt", root) ;
  if ((fp = fopenTag (root, "sites","r")))  { pbwtReadSites (p, fp) ; fclose (fp) ; }
  if ((fp = fopenTag (root, "samples","r"))) { pbwtReadSamples (p, fp) ; fclose (fp) ; }
  if ((fp = fopenTag (root, "missing","r"))) { pbwtReadMissing (p, fp) ; fclose (fp) ; }
  if ((fp = fopenTag (root, "dosage","r"))) { pbwtReadDosage (p, fp) ; fclose (fp) ; }
  if ((fp = fopenTag (root, "reverse","r"))) { pbwtReadReverse (p, fp) ; fclose (fp) ; }

  return p ;
}

/****************** MaCS file parsing *****************/

static void parseMacsHeader (FILE *fp, int *M, double *L)	/* parse MACS file header */
{
  if (strcmp (fgetword (fp), "COMMAND:")) die ("MaCS COMMAND line not found") ;
  fgetword (fp) ; 		/* the command */
  *M = atoi(fgetword(fp)) ; if (!*M) die ("failed to get M") ;
  *L = atof(fgetword(fp)) ; if (!*L) die ("failed to get L") ;
  while (fgetc(fp) != '\n') ;	/* ignore rest of line */
  
  if (strcmp (fgetword (fp), "SEED:")) die ("SEED line not found") ;
  while (fgetc(fp) != '\n') ;	/* ignore rest of line */
}

static BOOL parseMacsSite (FILE *fp, Site *s, int M, double L, uchar *yp) /* parse MACS site line */
{
  static uchar conv[256] ;
  static BOOL isFirst = TRUE ;
  int number ;

  if (isFirst) { conv['0'] = 0 ; conv['1'] = 1 ; isFirst = FALSE ; }
  if (feof (fp)) return FALSE ;

  if (strcmp (fgetword (fp), "SITE:")) return FALSE;

  number = atoi(fgetword(fp)) ;	/* this is the site number */
  s->x = (int) (L * atof(fgetword(fp))) ;
  atof(fgetword(fp)) ;		/* ignore the time */
  while (M--)
    *yp++ = conv[getc(fp)] ;
  if (feof (fp)) return FALSE ;
  if (getc(fp) != '\n') die ("end of line error for MaCS SITE %d", number) ;

  return TRUE ;
}

PBWT *pbwtReadMacs (FILE *fp)
{
  PBWT *p ;
  int M ;
  double L ;
  int i,j, nxPack, nyPack, nUnpack, n0 ;
  int nxTot = 0, nyTot = 0 ;
  uchar *x ;
  int *a ;
  PbwtCursor *u ;

  parseMacsHeader (fp, &M, &L) ;
  p = pbwtCreate (M, 0) ;
  p->sites = arrayCreate(4096, Site) ;
  u = pbwtCursorCreate (p, TRUE, TRUE) ;
  x = myalloc (M+1, uchar) ; x[M] = Y_SENTINEL ;	/* sentinel required for packing */

  while (parseMacsSite (fp, arrayp(p->sites,p->N,Site), M, L, x))
    { for (j = 0 ; j < M ; ++j) u->y[j] = x[u->a[j]] ; /* next character in sort order: BWT */
      pbwtCursorWriteForwards (u) ;
      p->N++ ;
      if (nCheckPoint && !(p->N % nCheckPoint)) pbwtCheckPoint (u, p) ;
    }
  pbwtCursorToAFend (u, p) ;

  fprintf (logFile, "read MaCS file: M, N are\t%d\t%d\n", M, p->N) ;
  if (isStats)
    fprintf (logFile, "                xtot, ytot are\t%d\t%d\n", nxTot, nyTot) ;

  free(x) ; pbwtCursorDestroy (u) ;

  return p ;
}

/************* read vcfq files, made with vcf query ... ***************/

static char* getVariation (FILE *fp)
{
  static Array a = 0 ;
  char c ;
  int i = 0 ;

  if (!a) a = arrayCreate (64, char) ;
  while (!isspace(c = getc(fp)) && !feof(fp)) array(a,i++,char) = c ;
  array(a,i++,char) = c ;
  while (!isspace(c = getc(fp)) && !feof(fp)) array(a,i++,char) = c ;
  array(a,i++,char) = 0 ;

  return arrp(a,0,char) ;
}

static BOOL parseVcfqLine (PBWT **pp, FILE *fp, Array x)
{
  char c, *chrom, *var ;
  int pos, m, len ;
  Site *s ;
  PBWT *p = *pp ;

  if (!p)			/* first call on a file - don't know M yet */
    chrom = strdup (fgetword (fp)) ;
  else
    while (!feof (fp) && !readMatchChrom (&p->chrom, fp))
      while (!feof (fp) && (getc (fp) != '\n')) ; /* ignore rest of line */
  if (feof (fp)) return 0 ;

  pos = atoi (fgetword (fp)) ;
  var = getVariation (fp) ;

  m = 0 ;
  while ((c = getc (fp)))
    switch (c)
      { case EOF: return 0 ;
      case '\n': goto endOfLine ; /* -1 is EOF */
      case '0': array(x,m++,uchar) = 0 ; break ;
      case '1': array(x,m++,uchar) = 1 ; break ;
      case '|': case '/': case '\\': case '\t': break ; /* could check ploidy here */
      default: die ("unexpected character %d in vcfq file genotype section") ;
      }
 endOfLine:
  if (feof (fp)) return FALSE ;
  if (p && m != p->M) die ("length mismatch reading vcfq line") ;

  if (!*pp)
    { p = *pp = pbwtCreate (m, 0) ;
      if (strcmp (chrom, ".")) p->chrom = chrom ;
      p->sites = arrayCreate(4096, Site) ;
      array(x,p->M,uchar) = Y_SENTINEL ; /* sentinel required for packing */
    }
  s = arrayp(p->sites, arrayMax(p->sites), Site) ;
  s->x = pos ;
  dictAdd (variationDict, var, &s->varD) ;

  ++p->N ;
  return TRUE ;
}

typedef BOOL (*ParseLineFunc)(PBWT** pp, FILE *f, Array a) ;

static PBWT *pbwtReadLineFile (FILE *fp, char* type, ParseLineFunc parseLine)
{
  PBWT *p = 0 ;
  int j ;
  uchar *x ;		/* original, sorted, compressed */
  int *a ;
  Array xArray = arrayCreate (10000, uchar) ;
  PbwtCursor *u ;

  while ((*parseLine) (&p, fp, xArray)) /* create p first time round */
    { if (!p->yz)		/* first line; p was just made! */
	{ p->yz = arrayCreate(4096*32, uchar) ;
	  u = pbwtCursorCreate (p, TRUE, TRUE) ;
	}
      x = arrp(xArray,0,uchar) ;
      for (j = 0 ; j < p->M ; ++j) u->y[j] = x[u->a[j]] ;
      pbwtCursorWriteForwards (u) ;
      if (nCheckPoint && !(p->N % nCheckPoint))	pbwtCheckPoint (u, p) ;
    }
  pbwtCursorToAFend (u, p) ;

  fprintf (logFile, "read %s file", type) ;
  if (p->chrom) fprintf (logFile, " for chromosome %s", p->chrom) ;
  fprintf (logFile, ": M, N are\t%d\t%d; yz length is %ld\n", p->M, p->N, arrayMax(p->yz)) ;

  arrayDestroy(xArray) ; pbwtCursorDestroy (u) ;

  return p ;
}

typedef BOOL (*ParseTwoLineFunc)(PBWT** pp, FILE *f, FILE *l, Array a) ;

static PBWT *pbwtReadLineTwoFile (FILE *fp, FILE *lp, char* type, ParseTwoLineFunc parseLine)
{
  PBWT *p = 0 ;
  int j ;
  uchar *x ;		/* original, sorted, compressed */
  int *a ;
  Array xArray = arrayCreate (10000, uchar) ;
  PbwtCursor *u ;

  /* skip header of legend file */
  while (!feof(lp))
  {  char c = getc(lp) ; if (c == '\n') break ;
     }

  while ((*parseLine) (&p, fp, lp, xArray)) /* create p first time round */
    { if (!p->yz)		/* first line; p was just made! */
	{ p->yz = arrayCreate(4096*32, uchar) ;
	  u = pbwtCursorCreate (p, TRUE, TRUE) ;
	}
      x = arrp(xArray,0,uchar) ;
      for (j = 0 ; j < p->M ; ++j) u->y[j] = x[u->a[j]] ;
      pbwtCursorWriteForwards (u) ;
      if (nCheckPoint && !(p->N % nCheckPoint))	pbwtCheckPoint (u, p) ;
    }
  pbwtCursorToAFend (u, p) ;

  fprintf (stderr, "read %s file", type) ;
  if (p->chrom) fprintf (stderr, " for chromosome %s", p->chrom) ;
  fprintf (stderr, ": M, N are\t%d\t%d; yz length is %ld\n", p->M, p->N, arrayMax(p->yz)) ;

  arrayDestroy(xArray) ; pbwtCursorDestroy (u) ;

  return p ;
}

PBWT *pbwtReadVcfq (FILE *fp) { return pbwtReadLineFile (fp, "vcfq", parseVcfqLine) ; }

/*************** read impute2 .gen format - contains sites not samples **********/

static long int nGenMissing ;

static BOOL parseGenLine (PBWT **pp, FILE *fp, Array x) /* based on parseVcfqLine() */
{
  PBWT *p = *pp ;

  fgetword(fp) ; fgetword(fp) ;	/* ignore first two name fields */

  int pos = atoi (fgetword (fp)) ;
  char *var = getVariation (fp) ; /* but need to change ' ' separator to '\t' */
  if (feof (fp)) return FALSE ;
  char *cp = var ; while (*cp && *cp != ' ') ++cp ; if (*cp == ' ') *cp = '\t' ; else die ("missing separator in line %d, var is %d", p?p->N:0, var) ;

  int m = 0, nscan ;
  while (!feof(fp))
    { char c = getc(fp) ; if (c == '\n') break ; else if (!isspace(c)) ungetc (c, fp) ;
      float f0, f1, f2 ;
      if ((nscan = fscanf (fp, "%f %f %f", &f0, &f1, &f2) != 3)) 
	die ("bad line %d, %d floats found, m %d, pos %d, var %s", p ? p->N : 0, nscan, m, pos, var) ;
      if (f0 + f1 + f2 == 0)	/* missing genotype */
	{ f0 = 1 ; ++nGenMissing ; }
      if (f0 + f1 + f2 < 0.98) 
	die ("inconsistent genotype in gen file: %f %f %f at %d line %d", 
	     f0, f1, f2, m, p ? p->N : 0) ;
      if (f0 > f1 && f0 > f2) { array(x,m++,uchar) = 0 ; array(x,m++,uchar) = 0 ; }
      else if (f1 > f2) { array(x,m++,uchar) = 0 ; array(x,m++,uchar) = 1 ; }
      else /* f2 is largest */ { array(x,m++,uchar) = 1 ; array(x,m++,uchar) = 1 ; }
    }

  if (feof (fp)) return FALSE ;
  if (p && m != p->M) die ("length mismatch reading vcfq line") ;

  if (!*pp)
    { p = *pp = pbwtCreate (m, 0) ;
      p->sites = arrayCreate(4096, Site) ;
      array(x,p->M,uchar) = Y_SENTINEL ; /* sentinel required for packing */
    }
  Site *s = arrayp(p->sites, arrayMax(p->sites), Site) ;
  s->x = pos ;
  dictAdd (variationDict, var, &s->varD) ;

  ++p->N ;
  return TRUE ;
}

static BOOL parseHapLine (PBWT **pp, FILE *fp, Array x) /* same as parseGenLine - slightly simpler */
{
  PBWT *p = *pp ;
  
  fgetword(fp) ;	fgetword(fp) ;	/* ignore first two name fields */
  
  int pos = atoi (fgetword (fp)) ;
  char *var = getVariation (fp) ; /* but need to change ' ' separator to '\t' */
  if (feof (fp)) return FALSE ;
  char *cp = var ; while (*cp && *cp != ' ') ++cp ; if (*cp == ' ') *cp = '\t' ; else die ("missing separator in line %d, var is %d", p?p->N:0, var) ;

  int m = 0, nscan ;
  while (!feof(fp))
  { char c = getc(fp) ; if (c == '\n') break ; else if (!isspace(c)) ungetc (c, fp) ;
    float f0, f1 ;
    if ((nscan = fscanf (fp, "%f %f", &f0, &f1) != 2))
      { warn ("bad line %d, %d floats found, m %d, pos %d, var %s - aborting", p ? p->N : 0, nscan, m, pos, var) ;
	return FALSE ;
      }
    array(x,m++,uchar) = f0 ; array(x,m++,uchar) = f1 ; /* haps are phased, put straight in as is */
  }
  
  if (feof (fp)) return FALSE ;
  if (p && m != p->M) die ("length mismatch reading haps line") ;
  
  if (!*pp)
  { p = *pp = pbwtCreate (m, 0) ;
    p->sites = arrayCreate(4096, Site) ;
    array(x,p->M,uchar) = Y_SENTINEL ; /* sentinel required for packing */
  }
  Site *s = arrayp(p->sites, arrayMax(p->sites), Site) ;
  s->x = pos ;
  dictAdd (variationDict, var, &s->varD) ;
  
  ++p->N ;
  return TRUE ;
}


static BOOL parseHapLegendLine (PBWT **pp, FILE *fp, FILE *lp, Array x) /* same as parseGenLine - slightly simpler */
{
  PBWT *p = *pp ;
  
  fgetword(lp) ;	/* ignore first name field */
  
  int pos = atoi (fgetword (lp)) ;
  char *var = getVariation (lp) ; /* but need to change ' ' separator to '\t' */
  if (feof (lp)) return FALSE ;
  char *cp = var ; while (*cp && *cp != ' ') ++cp ; if (*cp == ' ') *cp = '\t' ; else die ("missing separator in line %d, var is %d", p?p->N:0, var) ;
  while (!feof(lp))
  {  char c = getc(lp) ; if (c == '\n') break ;
     }

  int m = 0, nscan ;
  while (!feof(fp))
  { char c = getc(fp) ; if (c == '\n') break ; else if (!isspace(c)) ungetc (c, fp) ;
    float f0, f1 ;
    if ((nscan = fscanf (fp, "%f %f", &f0, &f1) != 2))
      { warn ("bad line %d, %d floats found, m %d, pos %d, var %s - aborting", p ? p->N : 0, nscan, m, pos, var) ;
	return FALSE ;
      }
    array(x,m++,uchar) = f0 ; array(x,m++,uchar) = f1 ; /* haps are phased, put straight in as is */
  }
  
  if (feof (fp)) return FALSE ;
  if (p && m != p->M) die ("length mismatch reading haps line") ;
  
  if (!*pp)
  { p = *pp = pbwtCreate (m, 0) ;
    p->sites = arrayCreate(4096, Site) ;
    array(x,p->M,uchar) = Y_SENTINEL ; /* sentinel required for packing */
  }
  Site *s = arrayp(p->sites, arrayMax(p->sites), Site) ;
  s->x = pos ;
  dictAdd (variationDict, var, &s->varD) ;
  
  ++p->N ;
  return TRUE ;
}

PBWT *pbwtReadGen (FILE *fp, char *chrom) 
{ 
  nGenMissing = 0 ;
  PBWT *p = pbwtReadLineFile (fp, "gen", parseGenLine) ;
  p->chrom = strdup (chrom) ;
  if (nGenMissing) fprintf (logFile, "%ld missing genotypes set to 00\n", nGenMissing) ;
  return p ;
}

PBWT *pbwtReadHap (FILE *fp, char *chrom)
{
  PBWT *p = pbwtReadLineFile (fp, "hap", parseHapLine) ;
  p->chrom = strdup (chrom) ;
  return p ;
}

PBWT *pbwtReadHapLegend (FILE *fp, FILE *lp, char *chrom)
{
  PBWT *p = pbwtReadLineTwoFile (fp, lp, "hap-legend", parseHapLegendLine) ;
  p->chrom = strdup (chrom) ;
  return p ;
}

PBWT *pbwtReadPhase (FILE *fp) /* Li and Stephens PHASE format */
{
  int nhaps=0,nsnps=0,ninds=0;
  int version=2;
  int l1=atoi(fgetword (fp)) ; if (getc(fp) != '\n') die ("bad first line in phase file") ;
  int l2 = atoi (fgetword (fp)) ;  if (getc(fp) != '\n') die ("bad second line in phase file") ;
  char * tc=fgetword (fp);
  int l3 = atoi (tc) ;  
  if (tc[0] == 'P') { // version 2
    nhaps=l1;
    nsnps=l2;
    ninds=nhaps/2;
  }else{ // version 1
    if (getc(fp) != '\n') die ("bad third line in phase file") ;
    ninds=l2;
    nhaps=l2*2;
    nsnps=l3;
    fgetword (fp); /* Remove the initial P */
    version=1;
  }
  fprintf(logFile,"Reading %i SNPs %i haplotypes and %i individuals from PHASE format version %i\n",nsnps,nhaps,ninds,version);
  PBWT *p = pbwtCreate (nhaps, nsnps) ;
  p->chrom = strdup ("0") ; /* chromosome number not available from phase files */
  p->sites = arrayCreate (4096, Site) ;
  int i ; for (i = 0 ; i < p->N ; ++i) {
    tc=fgetword(fp);
    arrayp(p->sites,i,Site)->x = atoi (tc) ;
  }
  if (getc(fp) != '\n') die ("bad location line in phase file") ;
  char var[2] ="S"; var[1] = 0 ; /*dummy variation line - used to create variationDict */
  for (i = 0 ; i < p->N ; ++i) 
    {
      if(version==1) *var = getc(fp) ;
      dictAdd (variationDict, var, &(arrayp(p->sites,i,Site)->varD)) ; }
  if(version==1) if (getc(fp) != '\n') die ("bad 5th line in phase file") ;

  uchar **data = myalloc (p->N, uchar*) ;
  for (i = 0 ; i < p->N ; ++i) data[i] = myalloc (p->M, uchar) ;
  int j ; for (j = 0 ; j < p->M ; ++j)
    { for (i = 0 ; i < p->N ; ++i) data[i][j] = getc(fp) - '0' ;
      if (getc(fp) != '\n') die ("bad %dth line in phase file", 7+j-version) ;
    }
  p->yz = arrayCreate(4096*32, uchar) ;
  PbwtCursor *u = pbwtCursorCreate (p, TRUE, TRUE) ;
  for (i = 0 ; i < p->N ; ++i)
    { for (j = 0 ; j < p->M ; ++j) u->y[j] = data[i][u->a[j]] ;
      pbwtCursorWriteForwards (u) ;
      if (nCheckPoint && !((i+1) % nCheckPoint)) pbwtCheckPoint (u, p) ;
    }
  pbwtCursorToAFend (u, p) ;

  fprintf (logFile, "read phase file") ;
  if (p->chrom) fprintf (logFile, " for chromosome %s", p->chrom) ;
  fprintf (logFile, ": M, N are\t%d\t%d; yz length is %ld\n", p->M, p->N, arrayMax(p->yz)) ;

  for (i = 0 ; i < p->N ; ++i) free(data[i]) ;
  free (data) ; pbwtCursorDestroy (u) ;
  return p ;
}

/*************** write haplotypes ******************/

void pbwtWriteHaplotypes (FILE *fp, PBWT *p)
{
  int i, j, n = 0, M = p->M ;
  uchar *hap = myalloc (M, uchar) ;
  PbwtCursor *u = pbwtCursorCreate (p, TRUE, TRUE) ;

  for (i = 0 ; i < p->N ; ++i)
    { for (j = 0 ; j < M ; ++j) hap[u->a[j]] = u->y[j] ;
      for (j = 0 ; j < M ; ++j) 
	{ if (isWriteImputeRef && j > 0) putc (' ', fp) ;
	  putc (hap[j]?'1':'0', fp) ;
	}
      putc ('\n', fp) ; fflush (fp) ;
      pbwtCursorForwardsRead (u) ;
    }
  free (hap) ; pbwtCursorDestroy (u) ;

  fprintf (logFile, "written haplotype file: %d rows of %d\n", p->N, M) ;
}

void pbwtWriteTransposedHaplotypes (PBWT *p, FILE *fp)
{
  int i, j ;
  uchar **hap = pbwtHaplotypes (p) ;
  
  for (j = 0 ; j < p->M ; ++j)
    { for (i = 0 ; i < p->N ; ++i) 
	putc (hap[j][i], fp) ;
      putc ('\n', fp) ; fflush (fp) ;
    }
  
  for (i = 0 ; j < p->M ; ++j) free(hap[j]) ;  free (hap) ;
  
  fprintf (logFile, "written transposed haplotype file: %d rows of %d\n", p->M,p->N) ;
}

/*************** write IMPUTE files ********************/

void pbwtWriteImputeRef (PBWT *p, char *fileNameRoot)
{
  FILE *fp ;

  isWriteImputeRef = TRUE ;

  if (!(fp=fopenTag (fileNameRoot, "imputeHaps","w"))) die ("can't open file") ; 
  pbwtWriteHaplotypes (fp, p) ; fclose (fp) ;

  if (!(fp=fopenTag (fileNameRoot, "imputeLegend","w"))) die ("can't open file") ; 
  fprintf (fp, "rsID\tposition\ta0\ta1\n") ; /* header line */
  pbwtWriteSites (p, fp) ; fclose (fp) ;

  isWriteImputeRef = FALSE ;
}

void pbwtWriteImputeHapsG (PBWT *p, FILE *fp)
{
  if (!p || !p->sites) die ("pbwtWriteImputeHaps called without sites") ;

  int i, j ;
  uchar *hap = myalloc (p->M, uchar) ;
  PbwtCursor *u = pbwtCursorCreate (p, TRUE, TRUE) ;

  for (i = 0 ; i < p->N ; ++i)
    { Site *s = arrp(p->sites, i, Site) ;
      fprintf (fp, "site%d\tsite%d\t%d", i+1, i+1, s->x) ;
      fprintf (fp, "\t%s", dictName (variationDict, s->varD)) ;
      
      for (j = 0 ; j < p->M ; ++j) hap[u->a[j]] = u->y[j] ;
      for (j = 0 ; j < p->M ; ++j) putc (' ', fp), putc (hap[j]?'1':'0', fp) ;
      putc ('\n', fp) ; fflush (fp) ;
      pbwtCursorForwardsRead (u) ;
    }

  free (hap) ; pbwtCursorDestroy (u) ;
}

void pbwtWriteGen (PBWT *p, FILE *fp)
{
  if (!p || !p->sites) die ("pbwtWriteImputeHaps called without sites") ;

  int i, j ;
  uchar *hap = myalloc (p->M, uchar) ;
  PbwtCursor *u = pbwtCursorCreate (p, TRUE, TRUE) ;
  BOOL isDosage = p->dosageOffset ? TRUE : FALSE ;
  double *d = 0, *ad = isDosage ? myalloc (p->M, double) : NULL;

  for (i = 0 ; i < p->N ; ++i)
    { Site *s = arrp(p->sites, i, Site) ;
      char *als = strdup( dictName(variationDict, s->varD) ), *ss = als ;
      while ( *ss ) { if ( *ss=='\t' ) *ss = '_' ; ss++ ; }
      fprintf (fp, "%s:%d_%s %s:%d_%s %d", p->chrom, s->x, als, p->chrom, s->x, als, s->x) ;
      ss = als ;
      while ( *ss ) { if ( *ss=='_' ) *ss = ' ' ; ss++ ; }
      fprintf (fp, " %s", als) ;
      free(als) ;
      if (isDosage) d = pbwtDosageRetrieve (p, u, d, i) ;
      
      for (j = 0 ; j < p->M ; ++j)
        { hap[u->a[j]] = u->y[j] ;
          if (isDosage) ad[u->a[j]] = d[j] ;
        }
      if (isDosage)
          for (j = 0 ; j < p->M ; j+=2) 
            fprintf (fp, " %f %f %f", (1-ad[j]) * (1-ad[j+1]), ad[j] + ad[j+1] - 2*ad[j]*ad[j+1], ad[j] * ad[j+1]) ;
      else
          for (j = 0 ; j < p->M ; j+=2) 
            if (hap[j] + hap[j+1] == 0) fprintf (fp, " 1 0 0") ;
            else if (hap[j] + hap[j+1] == 1) fprintf (fp, " 0 1 0") ;
            else fprintf (fp, " 0 0 1") ;
      putc ('\n', fp) ; fflush (fp) ;
      pbwtCursorForwardsRead (u) ;
    }

  if (ad) free(ad) ;
  free (hap) ; pbwtCursorDestroy (u) ;
}

/******************* end of file *******************/