text: /* psl.c was originally generated by the autoSql program, which also 
   * generated as_psl.h and as_psl.sql.  This module links the database and the RAM 
   * representation of objects. 
   *
   * This file is copyright 2002 Jim Kent, but license is hereby
   * granted for all use - public, private or commercial. */
  
  #include "common.h"
  #include "sqlNum.h"
  #include "sqlList.h"
  #include "localmem.h"
  #include "psl.h"
  #include "hash.h"
  #include "linefile.h"
  #include "dnaseq.h"
  #include "dystring.h"
  #include "fuzzyFind.h"
  #include "aliType.h"
  #include "binRange.h"
  #include "rangeTree.h"
  
  
  struct psl *pslxLoad(char **row)
  /* Load a psl from row fetched with select * from psl
   * from database.  Dispose of this with pslFree(). */
  {
  struct psl *ret = pslLoad(row);
  int retSize;
  sqlStringDynamicArray(row[21],&ret->qSequence, &retSize);
  sqlStringDynamicArray(row[22],&ret->tSequence, &retSize);
  return ret;
  }
  
  struct psl *pslLoad(char **row)
  /* Load a psl from row fetched with select * from psl
   * from database.  Dispose of this with pslFree(). */
  {
  struct psl *ret;
  int sizeOne;
  
  AllocVar(ret);
  ret->blockCount = sqlUnsigned(row[17]);
  ret->match = sqlUnsigned(row[0]);
  ret->misMatch = sqlUnsigned(row[1]);
  ret->repMatch = sqlUnsigned(row[2]);
  ret->nCount = sqlUnsigned(row[3]);
  ret->qNumInsert = sqlUnsigned(row[4]);
  ret->qBaseInsert = sqlSigned(row[5]);
  ret->tNumInsert = sqlUnsigned(row[6]);
  ret->tBaseInsert = sqlSigned(row[7]);
  strcpy(ret->strand, row[8]);
  ret->qName = cloneString(row[9]);
  ret->qSize = sqlUnsigned(row[10]);
  ret->qStart = sqlUnsigned(row[11]);
  ret->qEnd = sqlUnsigned(row[12]);
  ret->tName = cloneString(row[13]);
  ret->tSize = sqlUnsigned(row[14]);
  ret->tStart = sqlUnsigned(row[15]);
  ret->tEnd = sqlUnsigned(row[16]);
  sqlUnsignedDynamicArray(row[18], &ret->blockSizes, &sizeOne);
  if (sizeOne != ret->blockCount)
      {
      printf("sizeOne bloxksizes %d bs %d block=%s\n",sizeOne, ret->blockCount,row[18]);
      }
  assert(sizeOne == ret->blockCount);
  sqlUnsignedDynamicArray(row[19], &ret->qStarts, &sizeOne);
  if (sizeOne != ret->blockCount)
      {
      printf("sizeOne qStarts %d bs %d\n",sizeOne, ret->blockCount);
      }
  assert(sizeOne == ret->blockCount);
  sqlUnsignedDynamicArray(row[20], &ret->tStarts, &sizeOne);
  if (sizeOne != ret->blockCount)
      {
      printf("sizeOne tStarts %d bs %d\n",sizeOne, ret->blockCount);
      }
  assert(sizeOne == ret->blockCount);
  return ret;
  }
  
  struct psl *pslCommaIn(char **pS, struct psl *ret)
  /* Create a psl out of a comma separated string. 
   * This will fill in ret if non-null, otherwise will
   * return a new psl */
  {
  char *s = *pS;
  int i;
  
  if (ret == NULL)
      AllocVar(ret);
  ret->match = sqlUnsignedComma(&s);
  ret->misMatch = sqlUnsignedComma(&s);
  ret->repMatch = sqlUnsignedComma(&s);
  ret->nCount = sqlUnsignedComma(&s);
  ret->qNumInsert = sqlUnsignedComma(&s);
  ret->qBaseInsert = sqlSignedComma(&s);
  ret->tNumInsert = sqlUnsignedComma(&s);
  ret->tBaseInsert = sqlSignedComma(&s);
  sqlFixedStringComma(&s, ret->strand, sizeof(ret->strand));
  ret->qName = sqlStringComma(&s);
  ret->qSize = sqlUnsignedComma(&s);
  ret->qStart = sqlUnsignedComma(&s);
  ret->qEnd = sqlUnsignedComma(&s);
  ret->tName = sqlStringComma(&s);
  ret->tSize = sqlUnsignedComma(&s);
  ret->tStart = sqlUnsignedComma(&s);
  ret->tEnd = sqlUnsignedComma(&s);
  ret->blockCount = sqlUnsignedComma(&s);
  s = sqlEatChar(s, '{');
  AllocArray(ret->blockSizes, ret->blockCount);
  for (i=0; i<ret->blockCount; ++i)
      {
      ret->blockSizes[i] = sqlUnsignedComma(&s);
      }
  s = sqlEatChar(s, '}');
  s = sqlEatChar(s, ',');
  s = sqlEatChar(s, '{');
  AllocArray(ret->qStarts, ret->blockCount);
  for (i=0; i<ret->blockCount; ++i)
      {
      ret->qStarts[i] = sqlUnsignedComma(&s);
      }
  s = sqlEatChar(s, '}');
  s = sqlEatChar(s, ',');
  s = sqlEatChar(s, '{');
  AllocArray(ret->tStarts, ret->blockCount);
  for (i=0; i<ret->blockCount; ++i)
      {
      ret->tStarts[i] = sqlUnsignedComma(&s);
      }
  s = sqlEatChar(s, '}');
  s = sqlEatChar(s, ',');
  *pS = s;
  return ret;
  }
  
  void pslFree(struct psl **pEl)
  /* Free a single dynamically allocated psl such as created
   * with pslLoad(). */
  {
  struct psl *el;
  
  if ((el = *pEl) == NULL) return;
  freeMem(el->qName);
  freeMem(el->tName);
  freeMem(el->blockSizes);
  freeMem(el->qStarts);
  freeMem(el->tStarts);
  if (el->qSequence)
      {
      freeMem(el->qSequence[0]);
      freeMem(el->qSequence);
      }
  if (el->tSequence)
      {
      freeMem(el->tSequence[0]);
      freeMem(el->tSequence);
      }
  freez(pEl);
  }
  
  void pslFreeList(struct psl **pList)
  /* Free a list of dynamically allocated psl's */
  {
  struct psl *el, *next;
  
  for (el = *pList; el != NULL; el = next)
      {
      next = el->next;
      pslFree(&el);
      }
  *pList = NULL;
  }
  
  void pslOutputJson(struct psl *el, FILE *f) 
  /* print out psl as a json array */
  {
  fputs("[", f);
  fprintf(f, "%u,", el->match);
  fprintf(f, "%u,", el->misMatch);
  fprintf(f, "%u,", el->repMatch);
  fprintf(f, "%u,", el->nCount);
  fprintf(f, "%u,", el->qNumInsert);
  fprintf(f, "%d,", el->qBaseInsert);
  fprintf(f, "%u,", el->tNumInsert);
  fprintf(f, "%d,", el->tBaseInsert);
  fprintf(f, "\"%s\",", el->strand);
  fprintf(f, "\"%s\",", el->qName);
  fprintf(f, "%u,", el->qSize);
  fprintf(f, "%u,", el->qStart);
  fprintf(f, "%u,", el->qEnd);
  fprintf(f, "\"%s\",", el->tName);
  fprintf(f, "%u,", el->tSize);
  fprintf(f, "%u,", el->tStart);
  fprintf(f, "%u,", el->tEnd);
  fprintf(f, "%u,", el->blockCount);
  
  fputs("\"", f);
  for (int i=0; i<el->blockCount; ++i)
      {
      fprintf(f, "%u", el->blockSizes[i]);
      if (i < el->blockCount-1)
          fputc(',', f);
      }
  fputs("\",", f);
  
  fputs("\"", f);
  for (int i=0; i<el->blockCount; ++i)
      {
      fprintf(f, "%u", el->qStarts[i]);
      if (i < el->blockCount-1)
          fputc(',', f); // json does not allow trailing commas
      }
  fputs("\",", f);
  
  fputs("\"", f);
  for (int i=0; i<el->blockCount; ++i)
      {
      fprintf(f, "%u", el->tStarts[i]);
      if (i < el->blockCount-1)
          fputc(',', f);
      }
  fputs("\"", f);
  
  if (el->qSequence)
      {
      fputc(',',f);
      fputc('[',f);
      for (int i=0; i<el->blockCount; ++i)
  	{
  	fprintf(f, "'%s'", el->qSequence[i]);
          if (i < el->blockCount-1)
              fputc(',', f);
  	}
      fputc(']',f);
      fputc(',',f);
  
      fputc('[',f);
      for (int i=0; i<el->blockCount; ++i)
  	{
  	fprintf(f, "\"%s\"", el->tSequence[i]);
          if (i < el->blockCount-1)
              fputc(',', f);
  	}
      fputc(']',f);
      }
  
  if (ferror(f))
      {
      perror("Error writing psl file\n");
      errAbort("\n");
      }
  
  fputs("]", f);
  }
  
  void pslOutput(struct psl *el, FILE *f, char sep, char lastSep) 
  /* Print out psl.  Separate fields with sep. Follow last field with lastSep. */
  {
  int i;
  fprintf(f, "%u", el->match);
  fputc(sep,f);
  fprintf(f, "%u", el->misMatch);
  fputc(sep,f);
  fprintf(f, "%u", el->repMatch);
  fputc(sep,f);
  fprintf(f, "%u", el->nCount);
  fputc(sep,f);
  fprintf(f, "%u", el->qNumInsert);
  fputc(sep,f);
  fprintf(f, "%d", el->qBaseInsert);
  fputc(sep,f);
  fprintf(f, "%u", el->tNumInsert);
  fputc(sep,f);
  fprintf(f, "%d", el->tBaseInsert);
  fputc(sep,f);
  if (sep == ',') fputc('"',f);
  fprintf(f, "%s", el->strand);
  if (sep == ',') fputc('"',f);
  fputc(sep,f);
  if (sep == ',') fputc('"',f);
  fprintf(f, "%s", el->qName);
  if (sep == ',') fputc('"',f);
  fputc(sep,f);
  fprintf(f, "%u", el->qSize);
  fputc(sep,f);
  fprintf(f, "%u", el->qStart);
  fputc(sep,f);
  fprintf(f, "%u", el->qEnd);
  fputc(sep,f);
  if (sep == ',') fputc('"',f);
  fprintf(f, "%s", el->tName);
  if (sep == ',') fputc('"',f);
  fputc(sep,f);
  fprintf(f, "%u", el->tSize);
  fputc(sep,f);
  fprintf(f, "%u", el->tStart);
  fputc(sep,f);
  fprintf(f, "%u", el->tEnd);
  fputc(sep,f);
  fprintf(f, "%u", el->blockCount);
  fputc(sep,f);
  if (sep == ',') fputc('{',f);
  for (i=0; i<el->blockCount; ++i)
      {
      fprintf(f, "%u", el->blockSizes[i]);
      fputc(',', f);
      }
  if (sep == ',') fputc('}',f);
  fputc(sep,f);
  if (sep == ',') fputc('{',f);
  for (i=0; i<el->blockCount; ++i)
      {
      fprintf(f, "%u", el->qStarts[i]);
      fputc(',', f);
      }
  if (sep == ',') fputc('}',f);
  fputc(sep,f);
  if (sep == ',') fputc('{',f);
  for (i=0; i<el->blockCount; ++i)
      {
      fprintf(f, "%u", el->tStarts[i]);
      fputc(',', f);
      }
  if (sep == ',') fputc('}',f);
  if (el->qSequence)
      {
      fputc(sep,f);
      if (sep == ',') fputc('{',f);
      for (i=0; i<el->blockCount; ++i)
  	{
  	fprintf(f, "%s", el->qSequence[i]);
  	fputc(',', f);
  	}
      if (sep == ',') fputc('}',f);
      fputc(sep,f);
      if (sep == ',') fputc('{',f);
      for (i=0; i<el->blockCount; ++i)
  	{
  	fprintf(f, "%s", el->tSequence[i]);
  	fputc(',', f);
  	}
      if (sep == ',') fputc('}',f);
      }
  
  fputc(lastSep,f);
  if (ferror(f))
      {
      perror("Error writing psl file\n");
      errAbort("\n");
      }
  }
  
  /* ----- end autoSql generated part --------------- */
  
  void pslOutputShort(struct psl *el, FILE *f, char sep, char lastSep) 
  /* Print out psl.  Separate fields with sep. Follow last field with lastSep. */
  {
  fprintf(f, "%u", el->match);
  fputc(sep,f);
  fprintf(f, "%u", el->misMatch);
  fputc(sep,f);
  fprintf(f, "%u", el->repMatch);
  fputc(sep,f);
  fprintf(f, "%u", el->qNumInsert);
  fputc(sep,f);
  fprintf(f, "%d", el->qBaseInsert);
  fputc(sep,f);
  fprintf(f, "%u", el->tNumInsert);
  fputc(sep,f);
  fprintf(f, "%d", el->tBaseInsert);
  fputc(sep,f);
  if (sep == ',') fputc('"',f);
  fprintf(f, "%s", el->strand);
  if (sep == ',') fputc('"',f);
  fputc(sep,f);
  if (sep == ',') fputc('"',f);
  fprintf(f, "%s", el->qName);
  if (sep == ',') fputc('"',f);
  fputc(sep,f);
  fprintf(f, "%u", el->qStart);
  fputc(sep,f);
  fprintf(f, "%u", abs(el->qEnd - el->qStart));
  fputc(sep,f);
  if (sep == ',') fputc('"',f);
  fprintf(f, "%s", el->tName);
  if (sep == ',') fputc('"',f);
  fputc(sep,f);
  fprintf(f, "%u", el->tStart);
  fputc(sep,f);
  fprintf(f, "%u", abs(el->tEnd - el->tStart));
  fputc(sep,f);
  fprintf(f, "%u", el->blockCount);
  fputc(sep,f);
  if (sep == ',') fputc('{',f);
  fputc(lastSep,f);
  if (ferror(f))
      {
      perror("Error writing psl file\n");
      errAbort("\n");
      }
  }
  
  void pslOutFormat(struct psl *el, FILE *f, char sep, char lastSep) 
  /* Print out selected psl values.  Separate fields with sep. Follow last field with lastSep. */
  /* Prints out a better format with bold field headings followed by value */
  /* Requires further upstream work to ensure that only the field headers */
  /* declared here are printed if replacing an existing psl print function*/
  {
  const char *headers[] = {"Matches", "Mismatches", "Matches in repeats", "Number of N bases", "Query name", "Size", "Start", "End", "Chromosome", "Strand", "Start", "End"};
  char *hformat = "<B>%s:</B> "; /* string for formatted print for headers */
  char *uformat = "<B>%s:</B> %u%c"; /* string for formatted print for unsigned variable */
  char *targName;
  
  fprintf(f, uformat, headers[0], el->match, sep);
  fprintf(f, uformat, headers[1], el->misMatch, sep);
  fprintf(f, uformat, headers[2], el->repMatch, sep);
  fprintf(f, uformat, headers[3], el->nCount, sep);
  
  fprintf(f, hformat, headers[4]);
  if (sep == ',') fputc('"',f);
  fprintf(f, "%s", el->qName);
  if (sep == ',') fputc('"',f);
  fputc(sep,f);
  
  fprintf(f, uformat, headers[5], el->qSize, sep);
  fprintf(f, uformat, headers[6], el->qStart, sep);
  fprintf(f, uformat, headers[7], el->qEnd, sep);
  
  fprintf(f, hformat, headers[8]);
  if (sep == ',') fputc('"',f);
  /* skip leading 'chr' in string to get only chromosome part */
  targName = el->tName;
  if (startsWith("chr", el->tName))
     targName += 3;
  fprintf(f, "%s", targName);
  
  if (sep == ',') fputc('"',f);
  fputc(sep,f);
  
  fprintf(f, hformat, headers[9]);
  if (sep == ',') fputc('"',f);
  fprintf(f, "%s", el->strand);
  if (sep == ',') fputc('"',f);
  fputc(sep,f);
  
  fprintf(f, uformat, headers[10], el->tStart, sep);
  fprintf(f, uformat, headers[11], el->tEnd, sep);
  
  fputc(lastSep,f);
  
  if (ferror(f))
      {
      perror("Error writing psl file\n");
      errAbort("\n");
      }
  
  }
  
  struct psl *pslLoadAll(char *fileName)
  /* Load all psl's in file. */
  {
  struct lineFile *lf = pslFileOpen(fileName);
  struct psl *pslList = NULL, *psl;
  while ((psl = pslNext(lf)) != NULL)
      {
      slAddHead(&pslList, psl);
      }
  slReverse(&pslList);
  lineFileClose(&lf);
  return pslList;
  }
  
  
  int pslCmpQuery(const void *va, const void *vb)
  /* Compare to sort based on query start. */
  {
  const struct psl *a = *((struct psl **)va);
  const struct psl *b = *((struct psl **)vb);
  int dif;
  dif = strcmp(a->qName, b->qName);
  if (dif == 0)
      dif = a->qStart - b->qStart;
  return dif;
  }
  
  int pslCmpTarget(const void *va, const void *vb)
  /* Compare to sort based on target start. */
  {
  const struct psl *a = *((struct psl **)va);
  const struct psl *b = *((struct psl **)vb);
  int dif;
  dif = strcmp(a->tName, b->tName);
  if (dif == 0)
      dif = a->tStart - b->tStart;
  return dif;
  }
  
  int pslCmpTargetAndStrand(const void *va, const void *vb)
  /* Compare to sort based on target, strand,  tStart. */
  {
  const struct psl *a = *((struct psl **)va);
  const struct psl *b = *((struct psl **)vb);
  int dif;
  dif = strcmp(a->tName, b->tName);
  if (dif == 0)
      dif = strcmp(a->strand, b->strand);
  if (dif == 0)
      dif = a->tStart - b->tStart;
  return dif;
  }
  
  
  int pslCmpScore(const void *va, const void *vb)
  /* Compare to sort based on score (descending). */
  {
  const struct psl *a = *((struct psl **)va);
  const struct psl *b = *((struct psl **)vb);
  return pslScore(b) - pslScore(a);
  }
  
  int pslCmpQueryScore(const void *va, const void *vb)
  /* Compare to sort based on query then score (descending). */
  {
  const struct psl *a = *((struct psl **)va);
  const struct psl *b = *((struct psl **)vb);
  int diff = strcmp(a->qName, b->qName);
  if (diff == 0)
      diff = pslScore(b) - pslScore(a);
  return diff;
  }
  
  int pslCmpMatch(const void *va, const void *vb)
  /* Compare to sort based on match */
  {
  const struct psl *a = *((struct psl **)va);
  const struct psl *b = *((struct psl **)vb);
  return b->match - a->match;
  }
  
  static void pslLabelColumns(FILE *f)
  /* Write column info. */
  {
  fputs("\n"
  "match\tmis- \trep. \tN's\tQ gap\tQ gap\tT gap\tT gap\tstrand\tQ        \tQ   \tQ    \tQ  \tT        \tT   \tT    \tT  \tblock\tblockSizes \tqStarts\t tStarts\n"
  "     \tmatch\tmatch\t   \tcount\tbases\tcount\tbases\t      \tname     \tsize\tstart\tend\tname     \tsize\tstart\tend\tcount\n" 
  "---------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
  f);
  }
  
  static void pslLabelColumnsJson(FILE *f) 
  /* Write column info as a JSON array */
  {
  fputs("[\"matches\", \"misMatches\", \"repMatches\", \"nCount\", \"qNumInsert\", \"qBaseInsert\", "
          "\"tNumInsert\", \"tBaseInsert\", \"strand\", "
          "\"qName\", \"qSize\", \"qStart\", \"qEnd\", "
          "\"tName\", \"tSize\", \"tStart\", \"tEnd\", "
  	"\"blockCount\", \"blockSizes\", \"qStarts\", \"tStarts\"]", f);
  }
  
  void pslxWriteHead(FILE *f, enum gfType qType, enum gfType tType)
  /* Write header for extended (possibly protein) psl file. */
  {
  fprintf(f, "psLayout version 4 %s %s\n", gfTypeName(qType), gfTypeName(tType));
  pslLabelColumns(f);
  }
  
  void pslWriteHead(FILE *f)
  /* Write head of psl. */
  {
  fputs("psLayout version 3\n", f);
  pslLabelColumns(f);
  }
  
  void pslWriteAll(struct psl *pslList, char *fileName, boolean writeHeader)
  /* Write a psl file from list. */
  {
  FILE *f;
  struct psl *psl;
  
  f = mustOpen(fileName, "w");
  if (writeHeader)
      pslWriteHead(f);
  for (psl = pslList; psl != NULL; psl = psl->next)
      pslTabOut(psl, f);
  fclose(f);
  }
  
  void pslWriteAllJson(struct psl *pslList, FILE *f, char *db, boolean writeHeader)
  /* Write a psl file from list as a json array . */
  {
  fputs("{\n", f);
  if (writeHeader) 
      {
      fputs("\"track\": \"blat\",\n", f);
      fprintf(f, "\"genome\": \"%s\",\n", db);
      fputs("\"fields\": ", f);
      pslLabelColumnsJson(f);
      fputs(",\n", f);
      }
  fputs("\"blat\" : [\n", f);
  
  for (struct psl *psl = pslList; psl; psl = psl->next)
      {
      pslOutputJson(psl, f);
      if (psl->next)
          fputs(",\n", f);
      }
  
  puts("\n]\n}\n");
  }
  
  void pslxFileOpen(char *fileName, enum gfType *retQueryType, enum gfType *retTargetType, struct lineFile **retLf)
  /* Read header part of psl and make sure it's right.  Return
   * sequence types and file handle. */
  {
  char *line;
  int lineSize;
  char *words[30];
  char *version;
  int wordCount;
  int i;
  enum gfType qt = gftRna,  tt = gftDna;
  struct lineFile *lf = lineFileOpen(fileName, TRUE);
  
  if (!lineFileNext(lf, &line, &lineSize))
      warn("%s is empty", fileName);
  else
      {
      if (startsWith("psLayout version", line))
  	{
  	wordCount = chopLine(line, words);
  	if (wordCount < 3)
  	    errAbort("%s is not a psLayout file", fileName);
  	version = words[2];
  	if (sameString(version, "3"))
  	    {
  	    }
  	else if (sameString(version, "4"))
  	    {
  	    qt = gfTypeFromName(words[3]);
  	    tt = gfTypeFromName(words[4]);
  	    }
  	else
  	    {
  	    errAbort("%s is version %s of psLayout, this program can only handle through version 4",
  		fileName,  version);
  	    }
  	for (i=0; i<4; ++i)
  	    {
  	    if (!lineFileNext(lf, &line, &lineSize))
  		errAbort("%s severely truncated", fileName);
  	    }
  	}
      else
  	lineFileReuse(lf); 
      }
  *retQueryType = qt;
  *retTargetType = tt;
  *retLf = lf;
  }
  
  static void pslxFileOpenWithMetaConfig(char *fileName, bool isMetaUnique, enum gfType *retQueryType, enum gfType *retTargetType, struct lineFile **retLf, FILE *f)
  /* Read header part of psl and make sure it's right.  Return
   * sequence types and file handle and send meta data to output file f */
  {
  char *line;
  int lineSize;
  char *words[30];
  char *version;
  int wordCount;
  int i;
  enum gfType qt = gftRna,  tt = gftDna;
  struct lineFile *lf = lineFileOpen(fileName, TRUE);
  
  lineFileSetMetaDataOutput(lf, f);
  if (isMetaUnique)
      lineFileSetUniqueMetaData(lf);
  if (!lineFileNext(lf, &line, &lineSize))
      warn("%s is empty", fileName);
  else
      {
      if (startsWith("psLayout version", line))
  	{
  	wordCount = chopLine(line, words);
  	if (wordCount < 3)
  	    errAbort("%s is not a psLayout file", fileName);
  	version = words[2];
  	if (sameString(version, "3"))
  	    {
  	    }
  	else if (sameString(version, "4"))
  	    {
  	    qt = gfTypeFromName(words[3]);
  	    tt = gfTypeFromName(words[4]);
  	    }
  	else
  	    {
  	    errAbort("%s is version %s of psLayout, this program can only handle through version 4",
  		fileName,  version);
  	    }
  	for (i=0; i<4; ++i)
  	    {
  	    if (!lineFileNext(lf, &line, &lineSize))
  		errAbort("%s severely truncated", fileName);
  	    }
  	}
      else
          {
  	char *s = cloneString(line);
          boolean eof = FALSE;
          while ((line[0] == '#') && (!eof))
              {
              freeMem(s);
              if (!lineFileNext(lf, &line, &lineSize))
                  eof = TRUE;
              s = cloneString(line);
              }
  	wordCount = chopLine(s, words);
  	if ((wordCount < 21 || wordCount > 23 || (words[8][0] != '+' && words[8][0] != '-')) && (!eof))
  	    errAbort("%s is not a psLayout file", fileName);
  	else
  	    lineFileReuse(lf); 
  	freeMem(s);
  	}
      }
  *retQueryType = qt;
  *retTargetType = tt;
  *retLf = lf;
  }
  
  void pslxFileOpenWithMeta(char *fileName, enum gfType *retQueryType, enum gfType *retTargetType, struct lineFile **retLf, FILE *f)
  /* Read header part of psl and make sure it's right.  Return
   * sequence types and file handle and send meta data to output file f */
  {
  pslxFileOpenWithMetaConfig(fileName, FALSE, retQueryType, retTargetType, retLf, f);
  }
  
  void pslxFileOpenWithUniqueMeta(char *fileName, enum gfType *retQueryType, enum gfType *retTargetType, struct lineFile **retLf, FILE *f)
  /* Read header part of psl and make sure it's right.  Return
   * sequence types and file handle and send only unique meta data to output f */
  {
  pslxFileOpenWithMetaConfig(fileName, TRUE, retQueryType, retTargetType, retLf, f);
  }
  
  struct lineFile *pslFileOpen(char *fileName)
  /* Read header part of psl and make sure it's right. 
   * Return line file handle to it. */
  {
  enum gfType qt, tt;
  struct lineFile *lf;
  pslxFileOpen(fileName, &qt, &tt, &lf);
  return lf;
  }
  
  struct lineFile *pslFileOpenWithMeta(char *fileName, FILE *f)
  /* Read header part of psl and make sure it's right. 
   * Return line file handle to it. */
  {
  enum gfType qt, tt;
  struct lineFile *lf;
  pslxFileOpenWithMeta(fileName, &qt, &tt, &lf, f);
  return lf;
  }
  
  struct lineFile *pslFileOpenWithUniqueMeta(char *fileName, FILE *f)
  /* Read header part of psl and make sure it's right. 
   * Set flag to suppress duplicate header comments.
   * Return line file handle to it. */
  {
  enum gfType qt, tt;
  struct lineFile *lf;
  pslxFileOpenWithUniqueMeta(fileName, &qt, &tt, &lf, f);
  return lf;
  }
  
  struct psl *pslNext(struct lineFile *lf)
  /* Read next line from file and convert it to psl.  Return
   * NULL at eof. */
  {
  char *line;
  int lineSize;
  char *words[32];
  int wordCount;
  static int lineAlloc = 0;
  static char *chopBuf = NULL;
  
  if (!lineFileNextReal(lf, &line))
      {
      return NULL;
      }
  lineSize = strlen(line);
  if (lineSize >= lineAlloc)
      {
      lineAlloc = lineSize+256;
      chopBuf = needMoreMem(chopBuf, 0, lineAlloc);
      }
  memcpy(chopBuf, line, lineSize+1);
  wordCount = chopLine(chopBuf, words);
  if (wordCount == 21)
      {
      return pslLoad(words);
      }
  if (wordCount == 23)
      {
      return pslxLoad(words);
      }
  else
      {
      errAbort("Bad line %d of %s wordCount is %d instead of 21 or 23\n", lf->lineIx, lf->fileName, wordCount);
      return NULL;
      }
  }
  
  struct psl *pslxLoadLm(char **row, struct lm *lm)
  /* Load row into local memory pslx. */
  {
  struct psl *ret = pslLoadLm(row, lm);
  ret->qSequence = lmAlloc(lm, sizeof(ret->qSequence[0]) * ret->blockCount);
  sqlStringArray(lmCloneString(lm,row[21]),ret->qSequence, ret->blockCount);
  ret->tSequence = lmAlloc(lm, sizeof(ret->tSequence[0]) * ret->blockCount);
  sqlStringArray(lmCloneString(lm,row[22]),ret->tSequence, ret->blockCount);
  return ret;
  }
  
  struct psl *pslLoadLm(char **row, struct lm *lm)
  /* Load row into local memory psl. */
  {
  struct psl *ret;
  
  ret = lmAlloc(lm, sizeof(*ret));
  ret->blockCount = sqlUnsigned(row[17]);
  ret->match = sqlUnsigned(row[0]);
  ret->misMatch = sqlUnsigned(row[1]);
  ret->repMatch = sqlUnsigned(row[2]);
  ret->nCount = sqlUnsigned(row[3]);
  ret->qNumInsert = sqlUnsigned(row[4]);
  ret->qBaseInsert = sqlSigned(row[5]);
  ret->tNumInsert = sqlUnsigned(row[6]);
  ret->tBaseInsert = sqlSigned(row[7]);
  strcpy(ret->strand, row[8]);
  ret->qName = lmCloneString(lm,row[9]);
  ret->qSize = sqlUnsigned(row[10]);
  ret->qStart = sqlUnsigned(row[11]);
  ret->qEnd = sqlUnsigned(row[12]);
  ret->tName = lmCloneString(lm, row[13]);
  ret->tSize = sqlUnsigned(row[14]);
  ret->tStart = sqlUnsigned(row[15]);
  ret->tEnd = sqlUnsigned(row[16]);
  ret->blockSizes = lmAlloc(lm, sizeof(ret->blockSizes[0]) * ret->blockCount);
  sqlUnsignedArray(row[18], ret->blockSizes, ret->blockCount);
  ret->qStarts = lmAlloc(lm, sizeof(ret->qStarts[0]) * ret->blockCount);
  sqlUnsignedArray(row[19], ret->qStarts, ret->blockCount);
  ret->tStarts = lmAlloc(lm, sizeof(ret->tStarts[0]) * ret->blockCount);
  sqlUnsignedArray(row[20], ret->tStarts, ret->blockCount);
  return ret;
  }
  
  boolean pslIsProtein(const struct psl *psl)
  /* is psl a protein psl (are it's blockSizes and scores in protein space) */
  {
  int lastBlock = psl->blockCount - 1;
  
  return  (((psl->strand[1] == '+' ) &&
      (psl->tEnd == psl->tStarts[lastBlock] + 3*psl->blockSizes[lastBlock])) ||
     ((psl->strand[1] == '-') && 
      (psl->tStart == (psl->tSize-(psl->tStarts[lastBlock] + 3*psl->blockSizes[lastBlock])))));
  }
  
  int pslCalcMilliBad(struct psl *psl, boolean isMrna)
  /* Calculate badness in parts per thousand. */
  {
  int sizeMul = pslIsProtein(psl) ? 3 : 1;
  int qAliSize, tAliSize, aliSize;
  int milliBad = 0;
  int sizeDif;
  int insertFactor;
  int total;
  
  qAliSize = sizeMul * (psl->qEnd - psl->qStart);
  tAliSize = psl->tEnd - psl->tStart;
  aliSize = min(qAliSize, tAliSize);
  if (aliSize <= 0)
      return 0;
  sizeDif = qAliSize - tAliSize;
  if (sizeDif < 0)
      {
      if (isMrna)
  	sizeDif = 0;
      else
  	sizeDif = -sizeDif;
      }
  insertFactor = psl->qNumInsert;
  if (!isMrna)
      insertFactor += psl->tNumInsert;
  
  total = (sizeMul * (psl->match + psl->repMatch + psl->misMatch));
  if (total != 0)
      milliBad = (1000 * (psl->misMatch*sizeMul + insertFactor + round(3*log(1+sizeDif)))) / total;
  return milliBad;
  }
  
  int pslScore(const struct psl *psl)
  /* Return score for psl. */
  {
  int sizeMul = pslIsProtein(psl) ? 3 : 1;
  
  return sizeMul * (psl->match + ( psl->repMatch>>1)) - 
  	sizeMul * psl->misMatch - psl->qNumInsert - psl->tNumInsert;
  }
  
  int pslCmpScoreDesc(const void *va, const void *vb)
  /* Compare to sort based on score. */
  {
  const struct psl *a = *((struct psl **)va);
  const struct psl *b = *((struct psl **)vb);
  return pslScore(b) - pslScore(a);
  }
  
  
  struct ffAli *pslToFakeFfAli(struct psl *psl, DNA *needle, DNA *haystack)
  /* Convert from psl to ffAli format.  In some cases you can pass NULL
   * for needle and haystack - depending what the post-processing is going
   * to be. */
  {
  struct ffAli *ffList = NULL, *ff;
  int blockCount = psl->blockCount;
  unsigned *blockSizes = psl->blockSizes;
  unsigned *qStarts = psl->qStarts;
  unsigned *tStarts = psl->tStarts;
  int size;
  int i;
  
  for (i=0; i<blockCount; ++i)
      {
      size = blockSizes[i];
      AllocVar(ff);
      ff->left = ffList;
      ffList = ff;
      ff->nStart = ff->nEnd = needle + qStarts[i];
      ff->nEnd += size;
      ff->hStart = ff->hEnd = haystack + tStarts[i];
      ff->hEnd += size;
      }
  ffList = ffMakeRightLinks(ffList);
  return ffList;
  }
  
  struct psl *pslFromFakeFfAli(struct ffAli *ff, 
  	DNA *needle, DNA *haystack, char strand,
  	char *qName, int qSize, char *tName, int tSize)
  /* This will create a basic psl structure from a sorted series of ffAli
   * blocks.  The fields that would need actual sequence to be filled in
   * are left zero however - fields including match, repMatch, mismatch. */
  {
  struct psl *psl;
  unsigned *blockSizes;
  unsigned *qStarts;
  unsigned *tStarts;
  int blockCount;
  int i;
  int nStart, hStart;
  int nEnd, hEnd;
  
  AllocVar(psl);
  psl->blockCount = blockCount = ffAliCount(ff);
  psl->blockSizes = AllocArray(blockSizes, blockCount);
  psl->qStarts = AllocArray(qStarts, blockCount);
  psl->tStarts = AllocArray(tStarts, blockCount);
  psl->qName = cloneString(qName);
  psl->qSize = qSize;
  psl->tName = cloneString(tName);
  psl->tSize = tSize;
  psl->strand[0] = strand;
  
  for (i=0; i<blockCount; ++i)
      {
      nStart = ff->nStart - needle;
      nEnd = ff->nEnd - needle;
      hStart = ff->hStart - haystack;
      hEnd = ff->hEnd - haystack;
      blockSizes[i] = nEnd - nStart;
      qStarts[i] = nStart;
      tStarts[i] = hStart;
      if (i == 0)
         {
         psl->qStart = nStart;
         psl->tStart = hStart;
         }
      if (i == blockCount-1)
         {
         psl->qEnd = nEnd;
         psl->tEnd = hEnd;
         }
      ff = ff->right;
      }
  if (strand == '-')
      {
      reverseIntRange(&psl->qStart, &psl->qEnd, psl->qSize);
      }
  return psl;
  }
  
  struct ffAli *pslToFfAli(struct psl *psl, struct dnaSeq *query, struct dnaSeq *target,
  	int targetOffset)
  /* Convert from psl to ffAli format.  Clip to parts that we actually
   * have sequence for. */
  {
  struct ffAli *ffList = NULL, *ff;
  DNA *needle = query->dna;
  DNA *haystack = target->dna;
  int blockCount = psl->blockCount;
  unsigned *blockSizes = psl->blockSizes;
  unsigned *qStarts = psl->qStarts;
  unsigned *tStarts = psl->tStarts;
  int size;
  int i;
  int tMin = targetOffset;
  int tMax = targetOffset + target->size;
  int tStart, tEnd;
  int clipStart, clipEnd, clipOffset, clipSize;
  
  for (i=0; i<blockCount; ++i)
      {
      clipStart = tStart = tStarts[i];
      size = blockSizes[i];
      clipEnd = tEnd = tStart + size;
      if (tStart < tMax && tEnd > tMin)
  	{
  	if (clipStart < tMin) clipStart = tMin;
  	if (clipEnd > tMax) clipEnd = tMax;
  	clipOffset = clipStart - tStart;
  	clipSize = clipEnd - clipStart;
  	AllocVar(ff);
  	ff->left = ffList;
  	ffList = ff;
  	ff->nStart = ff->nEnd = needle + qStarts[i] + clipOffset;
  	ff->nEnd += clipSize;
  	ff->hStart = ff->hEnd = haystack + clipStart - targetOffset;
  	ff->hEnd += clipSize;
  	}
      }
  ffList = ffMakeRightLinks(ffList);
  ffCountGoodEnds(ffList);
  return ffList;
  }
  
  int pslOrientation(struct psl *psl)
  /* Translate psl strand + or - to orientation +1 or -1 */
  {
  if (psl->strand[1] != '\0')
      {
      /* translated blat */
      if (psl->strand[0] != psl->strand[1])
          return -1;
      else
          return 1;
      }
  else
      {
      if (psl->strand[0] == '-')
          return -1;
      else
          return 1;
      }
  }
  
  int pslWeightedIntronOrientation(struct psl *psl, struct dnaSeq *genoSeq, int offset)
  /* Return >0 if introns make it look like alignment is on + strand,
   *        <0 if introns make it look like alignment is on - strand,
   *        0 if can't tell.  The absolute value of the return indicates
   * how many splice sites we've seen supporting the orientation.
   * Sequence should NOT be reverse complemented.  */
  {
  int intronDir = 0;
  int oneDir;
  int i;
  DNA *dna = genoSeq->dna;
  
  /* code below doesn't support negative target strand (translated blat) */
  if (psl->strand[1] == '-')
      errAbort("pslWeightedIntronOrientation doesn't support a negative target strand");
  
  for (i=1; i<psl->blockCount; ++i)
      {
      int iStart, iEnd, blockSize = psl->blockSizes[i-1];
      if (psl->qStarts[i-1] + blockSize == psl->qStarts[i])
  	{
  	iStart = psl->tStarts[i-1] + psl->blockSizes[i-1] - offset;
  	iEnd = psl->tStarts[i] - offset;
  	oneDir = intronOrientation(dna+iStart, dna+iEnd);
  	intronDir += oneDir;
  	}
      }
  return intronDir;
  }
  
  int pslIntronOrientation(struct psl *psl, struct dnaSeq *genoSeq, int offset)
  /* Return 1 if introns make it look like alignment is on + strand,
   *       -1 if introns make it look like alignment is on - strand,
   *        0 if can't tell.
   * Sequence should NOT be reverse complemented.  */
  {
  int intronDir = pslWeightedIntronOrientation(psl, genoSeq, offset);
  if (intronDir < 0)
      intronDir = -1;
  else if (intronDir > 0)
      intronDir = 1;
  return intronDir;
  }
  
  boolean pslHasIntron(struct psl *psl, struct dnaSeq *seq, int seqOffset)
  /* Return TRUE if there's a probable intron. Sequence should NOT be
   * reverse complemented.*/
  {
  int blockCount = psl->blockCount, i;
  unsigned *tStarts = psl->tStarts;
  unsigned *blockSizes = psl->blockSizes;
  unsigned *qStarts = psl->qStarts;
  int blockSize, start, end;
  DNA *dna = seq->dna;
  
  for (i=1; i<blockCount; ++i)
      {
      blockSize = blockSizes[i-1];
      start = qStarts[i-1]+blockSize;
      end = qStarts[i];
      if (start == end)
          {
          start = tStarts[i-1] + blockSize;
          end = tStarts[i];
          if (psl->strand[1] == '-')
              reverseIntRange(&start, &end, psl->tSize);
          start -= seqOffset;
          end -= seqOffset;
  	if (intronOrientation(dna+start, dna+end) != 0)
  	    return TRUE;
  	}
      }
  return FALSE;
  }
  
  void pslTailSizes(struct psl *psl, int *retStartTail, int *retEndTail)
  /* Find the length of "tails" (rather than extensions) implied by psl. */
  {
  int orientation = pslOrientation(psl);
  int qFloppyStart, qFloppyEnd;
  int tFloppyStart, tFloppyEnd;
  
  if (orientation > 0)
      {
      qFloppyStart = psl->qStart;
      qFloppyEnd = psl->qSize - psl->qEnd;
      }
  else
      {
      qFloppyStart = psl->qSize - psl->qEnd;
      qFloppyEnd = psl->qStart;
      }
  tFloppyStart = psl->tStart;
  tFloppyEnd = psl->tSize - psl->tEnd;
  *retStartTail = min(qFloppyStart, tFloppyStart);
  *retEndTail = min(qFloppyEnd, tFloppyEnd);
  }
  
  static void rcSeqs(char **seqs, unsigned blockCount, unsigned *blockSizes)
  /* reverses complement sequences in list, maintain property that all strings
   * are in one malloc block.   blockSizes should already be reversed. */
  {
  char *buf, *next;
  int i, memSz = 0;
  
  /* get a new memory block for strings */
  for (i = 0; i < blockCount; i++)
      memSz += blockSizes[i]+1;
  next = buf = needLargeMem(memSz);
  
  /* reverse compliment and copy to new memory block */
  for (i = blockCount-1; i >= 0; i--)
      {
      int len = strlen(seqs[i]);
      reverseComplement(seqs[i], len);
      memcpy(next, seqs[i], len+1);
      next += len+1;
      }
  
  /* swap memory and update pointers */
  freeMem(seqs[0]);
  seqs[0] = buf;
  next = buf;
  
  for (i = 0; i < blockCount; i++)
      {
      seqs[i] = next;
      next += blockSizes[i]+1;
      }
  }
  
  void pslRc(struct psl *psl)
  /* Reverse-complement a PSL alignment.  This makes the target strand explicit. */
  {
  unsigned tSize = psl->tSize, qSize = psl->qSize;
  unsigned blockCount = psl->blockCount, i;
  unsigned *tStarts = psl->tStarts, *qStarts = psl->qStarts, *blockSizes = psl->blockSizes;
  int mult = pslIsProtein(psl) ? 3 : 1;
  
  /* swap strand, forcing target to have an explict strand */
  psl->strand[0] = (psl->strand[0] != '-') ? '-' : '+';
  psl->strand[1] = (psl->strand[1] != '-') ? '-' : '+';
  psl->strand[2] = 0;
  
  for (i=0; i<blockCount; ++i)
      {
      tStarts[i] = tSize - (tStarts[i] + mult * blockSizes[i]);
      qStarts[i] = qSize - (qStarts[i] + blockSizes[i]);
      }
  reverseUnsigned(tStarts, blockCount);
  reverseUnsigned(qStarts, blockCount);
  reverseUnsigned(blockSizes, blockCount);
  if (psl->qSequence != NULL)
      {
      rcSeqs(psl->qSequence, blockCount, blockSizes);
      rcSeqs(psl->tSequence, blockCount, blockSizes);
      }
  }
  
  
  /* macro to swap to variables */
  #define swapVars(a, b, tmp) ((tmp) = (a), (a) = (b), (b) = (tmp))
  
  static void swapBlocks(struct psl *psl)
  /* Swap the blocks in a psl without reverse complementing them. */
  {
  int i;
  unsigned utmp;
  char *stmp; 
  for (i = 0; i < psl->blockCount; i++)
      {
      swapVars(psl->qStarts[i], psl->tStarts[i], utmp);
      if (psl->qSequence != NULL)
          swapVars(psl->qSequence[i], psl->tSequence[i], stmp);
      }
  }
  
  static void swapRCBlocks(struct psl *psl)
  /* Swap and reverse complement blocks in a psl. Other psl fields must
   * be modified first */
  {
  int i;
  unsigned *uatmp;
  char **satmp;
  reverseUnsigned(psl->tStarts, psl->blockCount);
  reverseUnsigned(psl->qStarts, psl->blockCount);
  reverseUnsigned(psl->blockSizes, psl->blockCount);
  swapVars(psl->tStarts, psl->qStarts, uatmp);
  
  /* qSize and tSize have already been swapped */
  for (i = 0; i < psl->blockCount; i++)
      {
      psl->qStarts[i] = psl->qSize - (psl->qStarts[i] + psl->blockSizes[i]);
      psl->tStarts[i] = psl->tSize - (psl->tStarts[i] + psl->blockSizes[i]);
      }
  if (psl->qSequence != NULL)
      {
      /* note: all block sequences are stored in one malloc block, which is
       * entry zero */
      rcSeqs(psl->qSequence, psl->blockCount, psl->blockSizes);
      rcSeqs(psl->tSequence, psl->blockCount, psl->blockSizes);
      swapVars(psl->qSequence, psl->tSequence, satmp);
      }
  }
  
  void pslSwap(struct psl *psl, boolean noRc)
  /* swap query and target in psl.  If noRc is TRUE, don't reverse-complement
   * PSL if needed, instead make target strand explict. */
  {
  int itmp;
  unsigned utmp;
  char ctmp, *stmp; 
  swapVars(psl->qBaseInsert, psl->tBaseInsert, utmp);
  swapVars(psl->tNumInsert, psl->qNumInsert, utmp);
  swapVars(psl->qName, psl->tName, stmp);
  swapVars(psl->qSize, psl->tSize, utmp);
  swapVars(psl->qStart, psl->tStart, itmp);
  swapVars(psl->qEnd, psl->tEnd, itmp);
  
  /* handle strand and block copy */
  if (psl->strand[1] != '\0')
      {
      /* translated */
      swapVars(psl->strand[0], psl->strand[1], ctmp);
      swapBlocks(psl);
      }
  else if (noRc)
      {
      /* untranslated with no reverse complement */
      psl->strand[1] = psl->strand[0];
      psl->strand[0] = '+';
      swapBlocks(psl);
      }
  else
      {
      /* untranslated */
      if (psl->strand[0] == '+')
          swapBlocks(psl);
      else
          swapRCBlocks(psl);
      }
  }
  
  void pslTargetOffset(struct psl *psl, int offset)
  /* Add offset to target positions in psl. */
  {
  int i, blockCount = psl->blockCount;
  unsigned *tStarts = psl->tStarts;
  psl->tStart += offset;
  psl->tEnd += offset;
  for (i=0; i<blockCount; ++i)
     tStarts[i] += offset;
  }
  
  void pslDump(struct psl *psl, FILE *f)
  /* Dump most of PSL to file - for debugging. */
  {
  int i;
  fprintf(f, "<PRE>\n");
  fprintf(f, "psl %s:%d-%d %s %s:%d-%d %d\n", 
  	psl->qName, psl->qStart, psl->qEnd, psl->strand,
  	psl->tName, psl->tStart, psl->tEnd, psl->blockCount);
  for (i=0; i<psl->blockCount; ++i)
      fprintf(f, "  size %d, qStart %d, tStart %d\n", 
      	psl->blockSizes[i], psl->qStarts[i], psl->tStarts[i]);
  fprintf(f, "</PRE>");
  }
  
  void pslRecalcBounds(struct psl *psl)
  /* Calculate qStart/qEnd tStart/tEnd at top level to be consistent
   * with blocks. */
  {
  int qStart, qEnd, tStart, tEnd, size;
  int last = psl->blockCount - 1;
  qStart = psl->qStarts[0];
  tStart = psl->tStarts[0];
  size = psl->blockSizes[last];
  qEnd = psl->qStarts[last] + size;
  tEnd = psl->tStarts[last] + size;
  if (psl->strand[0] == '-')
      reverseIntRange(&qStart, &qEnd, psl->qSize);
  if (psl->strand[1] == '-')
      reverseIntRange(&tStart, &tEnd, psl->tSize);
  psl->qStart = qStart;
  psl->qEnd = qEnd;
  psl->tStart = tStart;
  psl->tEnd = tEnd;
  }
  
  struct psl *pslTrimToTargetRange(struct psl *oldPsl, int tMin, int tMax)
  /* Return psl trimmed to fit inside tMin/tMax.  Note this does not
   * update the match/misMatch and related fields. */
  {
  int newSize;
  int oldBlockCount = oldPsl->blockCount;
  boolean tIsRc = (oldPsl->strand[1] == '-');
  int newBlockCount = 0, completeBlockCount = 0;
  int i;
  struct psl *newPsl = NULL;
  int tMn = tMin, tMx = tMax;   /* tMin/tMax adjusted for strand. */
  
  /* Deal with case where we're completely trimmed out quickly. */
  newSize = rangeIntersection(oldPsl->tStart, oldPsl->tEnd, tMin, tMax);
  if (newSize <= 0)
      return NULL;
  
  if (tIsRc)
      reverseIntRange(&tMn, &tMx, oldPsl->tSize);
  
  /* Count how many blocks will survive trimming. */
  oldBlockCount = oldPsl->blockCount;
  for (i=0; i<oldBlockCount; ++i)
      {
      int s = oldPsl->tStarts[i];
      int e = s + oldPsl->blockSizes[i];
      int sz = e - s;
      int overlap;
      if ((overlap = rangeIntersection(s, e, tMn, tMx)) > 0)
          ++newBlockCount;
      if (overlap == sz)
          ++completeBlockCount;
      }
  
  if (newBlockCount == 0)
      return NULL;
  
  /* Allocate new psl and fill in what we already know. */
  AllocVar(newPsl);
  strcpy(newPsl->strand, oldPsl->strand);
  newPsl->qName = cloneString(oldPsl->qName);
  newPsl->qSize = oldPsl->qSize;
  newPsl->tName = cloneString(oldPsl->tName);
  newPsl->tSize = oldPsl->tSize;
  newPsl->blockCount = newBlockCount;
  AllocArray(newPsl->blockSizes, newBlockCount);
  AllocArray(newPsl->qStarts, newBlockCount);
  AllocArray(newPsl->tStarts, newBlockCount);
  
  /* Fill in blockSizes, qStarts, tStarts with real data. */
  newBlockCount = completeBlockCount = 0;
  for (i=0; i<oldBlockCount; ++i)
      {
      int oldSz = oldPsl->blockSizes[i];
      int sz = oldSz;
      int tS = oldPsl->tStarts[i];
      int tE = tS + sz;
      int qS = oldPsl->qStarts[i];
      int qE = qS + sz;
      if (rangeIntersection(tS, tE, tMn, tMx) > 0)
          {
  	int diff;
  	if ((diff = (tMn - tS)) > 0)
  	    {
  	    tS += diff;
  	    qS += diff;
  	    sz -= diff;
  	    }
  	if ((diff = (tE - tMx)) > 0)
  	    {
  	    tE -= diff;
  	    qE -= diff;
  	    sz -= diff;
  	    }
  	newPsl->qStarts[newBlockCount] = qS;
  	newPsl->tStarts[newBlockCount] = tS;
  	newPsl->blockSizes[newBlockCount] = sz;
  	++newBlockCount;
  	if (sz == oldSz)
  	    ++completeBlockCount;
  	}
      }
  pslRecalcBounds(newPsl);
  return newPsl;
  }
  
  struct psl *pslTrimToQueryRange(struct psl *oldPsl, int qMin, int qMax)
  /* Return psl trimmed to fit inside qMin/qMax.  Note this does not
   * update the match/misMatch and related fields. */
  {
  int newSize;
  int oldBlockCount = oldPsl->blockCount;
  boolean qIsRc = (oldPsl->strand[0] == '-');
  int newBlockCount = 0, completeBlockCount = 0;
  int i;
  struct psl *newPsl = NULL;
  int qMn = qMin, qMx = qMax;   /* qMin/qMax adjusted for strand. */
  
  /* Deal with case where we're completely trimmed out quickly. */
  newSize = rangeIntersection(oldPsl->qStart, oldPsl->qEnd, qMin, qMax);
  if (newSize <= 0)
      return NULL;
  
  if (qIsRc)
      reverseIntRange(&qMn, &qMx, oldPsl->qSize);
  
  /* Count how many blocks will survive trimming. */
  oldBlockCount = oldPsl->blockCount;
  for (i=0; i<oldBlockCount; ++i)
      {
      int s = oldPsl->qStarts[i];
      int e = s + oldPsl->blockSizes[i];
      int sz = e - s;
      int overlap;
      if ((overlap = rangeIntersection(s, e, qMn, qMx)) > 0)
          ++newBlockCount;
      if (overlap == sz)
          ++completeBlockCount;
      }
  
  if (newBlockCount == 0)
      return NULL;
  
  /* Allocate new psl and fill in what we already know. */
  AllocVar(newPsl);
  strcpy(newPsl->strand, oldPsl->strand);
  newPsl->qName = cloneString(oldPsl->qName);
  newPsl->qSize = oldPsl->qSize;
  newPsl->tName = cloneString(oldPsl->tName);
  newPsl->tSize = oldPsl->tSize;
  newPsl->blockCount = newBlockCount;
  AllocArray(newPsl->blockSizes, newBlockCount);
  AllocArray(newPsl->qStarts, newBlockCount);
  AllocArray(newPsl->tStarts, newBlockCount);
  
  /* Fill in blockSizes, qStarts, tStarts with real data. */
  newBlockCount = completeBlockCount = 0;
  for (i=0; i<oldBlockCount; ++i)
      {
      int oldSz = oldPsl->blockSizes[i];
      int sz = oldSz;
      int qS = oldPsl->qStarts[i];
      int qE = qS + sz;
      int tS = oldPsl->tStarts[i];
      int tE = tS + sz;
      if (rangeIntersection(qS, qE, qMn, qMx) > 0)
          {
  	int diff;
  	if ((diff = (qMn - qS)) > 0)
  	    {
  	    tS += diff;
  	    qS += diff;
  	    sz -= diff;
  	    }
  	if ((diff = (qE - qMx)) > 0)
  	    {
  	    tE -= diff;
  	    qE -= diff;
  	    sz -= diff;
  	    }
  	newPsl->qStarts[newBlockCount] = qS;
  	newPsl->tStarts[newBlockCount] = tS;
  	newPsl->blockSizes[newBlockCount] = sz;
  	++newBlockCount;
  	if (sz == oldSz)
  	    ++completeBlockCount;
  	}
      }
  pslRecalcBounds(newPsl);
  return newPsl;
  }
  
  static void printPslDesc(char* pslDesc, FILE* out, struct psl* psl)
  /* print description of a PSL on first error */
  {
  fprintf(out, "Error: invalid PSL: %s:%u-%u %s:%u-%u %s %s\n",
          psl->qName, psl->qStart, psl->qEnd,
          psl->tName, psl->tStart, psl->tEnd,
          psl->strand, pslDesc);
  }
  
  
  static void chkError(char* pslDesc, FILE* out, struct psl* psl, int* errCount, char* format, ...)
  /* forward needed to specify printf signature for gcc checking */
  #if defined(__GNUC__)
  __attribute__((format(printf, 5, 6)))
  #endif
  ;
  
  static void chkError(char* pslDesc, FILE* out, struct psl* psl, int* errCount, char* format, ...)
  /* error handling on an pslCheck error, counting error and issuing description
   * of PSL on the first error. */
  {
  if (*errCount == 0)
      printPslDesc(pslDesc, out, psl);
  va_list args;
  va_start(args, format);
  vfprintf(out, format, args);
  va_end(args);
  (*errCount)++;
  }
  
  static void chkBlkRanges(char* pslDesc, FILE* out, struct psl* psl,
                           char* pName, char* pLabel, char pCLabel, char pStrand,
                           unsigned pSize, unsigned pStart, unsigned pEnd,
                           unsigned iBlk, unsigned* pBlockStarts, int* errCount)
  /* check the target or query ranges in a PSL incrementing errorCnt */
  {
  unsigned blkStart = pBlockStarts[iBlk];
  unsigned blkEnd = blkStart+psl->blockSizes[iBlk];
  /* translate stand to genomic coords */
  unsigned gBlkStart = (pStrand == '+') ? blkStart : (pSize - blkEnd);
  unsigned gBlkEnd = (pStrand == '+') ? blkEnd : (pSize - blkStart);
  
  if ((pSize > 0) && (blkEnd > pSize))
      chkError(pslDesc, out, psl, errCount,
               "\t%s %s block %u end %u > %cSize %u\n",
               pName, pLabel, iBlk, blkEnd, pCLabel, pSize);
  if (gBlkStart < pStart)
      chkError(pslDesc, out, psl, errCount,
               "\t%s %s block %u start %u < %cStart %u\n",
               pName, pLabel, iBlk, gBlkStart, pCLabel, pStart);
  if (gBlkStart >= pEnd)
      chkError(pslDesc, out, psl, errCount,
               "\t%s %s block %u start %u >= %cEnd %u\n",
               pName, pLabel, iBlk, gBlkStart, pCLabel, pEnd);
  if (gBlkEnd < pStart)
      chkError(pslDesc, out, psl, errCount,
               "\t%s %s block %u end %u < %cStart %u\n",
               pName, pLabel, iBlk, gBlkEnd, pCLabel, pStart);
  if (gBlkEnd > pEnd)
      chkError(pslDesc, out, psl, errCount,
               "\t%s %s block %u end %u > %cEnd %u\n",
               pName, pLabel, iBlk, gBlkEnd, pCLabel, pEnd);
  if (iBlk > 0)
      {
      unsigned prevBlkEnd = pBlockStarts[iBlk-1]+psl->blockSizes[iBlk-1];
      if (blkStart < prevBlkEnd)
          chkError(pslDesc, out, psl, errCount,
                   "\t%s %s block %u start %u < previous block end %u\n",
                   pName, pLabel, iBlk, blkStart, prevBlkEnd);
      }
  }
  
  static void chkRanges(char* pslDesc, FILE* out, struct psl* psl,
                        char* pName, char* pLabel, char pCLabel, char pStrand,
                        unsigned pSize, unsigned pStart, unsigned pEnd,
                        unsigned* pBlockStarts, int blockSizeMult, int* errCount)
  /* check the target or query ranges in a PSL, increment errorCnt */
  {
  unsigned iBlk;
  if (pStart >= pEnd)
      chkError(pslDesc, out, psl, errCount,
               "\t%s %cStart %u >= %cEnd %u\n",
               pName, pCLabel, pStart, pCLabel, pEnd);
  if (pEnd > pSize)
      chkError(pslDesc, out, psl, errCount,
               "\t%s %cEnd %u >= %cSize %u\n",
               pName, pCLabel, pEnd, pCLabel, pSize);
  // check that block start/end matches overall start end
  unsigned pStartStrand = pStart, pEndStrand = pEnd;
  if (pStrand != '+')
      reverseUnsignedRange(&pStartStrand, &pEndStrand, pSize);
  unsigned lastBlkEnd = pBlockStarts[psl->blockCount-1] + (blockSizeMult * psl->blockSizes[psl->blockCount-1]);
  if ((pStartStrand != pBlockStarts[0]) || (pEndStrand != lastBlkEnd))
      chkError(pslDesc, out, psl, errCount,
               "\t%s strand \"%c\" adjusted %cStart-%cEnd range %u-%u != block range %u-%u\n",
               pName, pStrand, pCLabel, pCLabel, pStartStrand, pEndStrand, pBlockStarts[0], lastBlkEnd);
  
  for (iBlk = 0; iBlk < psl->blockCount; iBlk++)
      chkBlkRanges(pslDesc, out, psl, pName, pLabel, pCLabel, pStrand,
                   pSize, pStart, pEnd, iBlk, pBlockStarts, errCount);
  }
  
  
  static void chkInsertCounts(char* pslDesc, FILE* out, struct psl* psl,
                              char* pName, char pCLabel, unsigned* pBlockStarts,
                              unsigned pNumInsert, unsigned pBaseInsert,
                              int* errCount)
  /* check the insert counts, incrementing errorCnt */
  {
  unsigned numInsert = 0, baseInsert = 0;
  int iBlk;
  
  for (iBlk = 1; iBlk < psl->blockCount; iBlk++)
      {
      unsigned gapSize = pBlockStarts[iBlk] - (pBlockStarts[iBlk-1]+psl->blockSizes[iBlk-1]);
      if (gapSize > 0)
          {
          numInsert++;
          baseInsert += gapSize;
          }
      }
  if (numInsert != pNumInsert)
      chkError(pslDesc, out, psl, errCount,
               "\t%s %cNumInsert %u != expected %u\n",
               pName, pCLabel, pNumInsert, numInsert);
  if (baseInsert != pBaseInsert)
      chkError(pslDesc, out, psl, errCount,
               "\t%s %cBaseInsert %u != expected %u\n",
               pName, pCLabel, pBaseInsert, baseInsert);
  }
  
  int pslCheck2(unsigned opts, char *pslDesc, FILE* out, struct psl* psl)
  /* Validate a PSL for consistency.  pslDesc is printed the error messages to
   * file out (open /dev/null to discard). Return count of errors.  Option
   * PSL_CHECK_IGNORE_INSERT_CNTS doesn't validate problems insert counts fields
   * in each PSL.  Useful because protein PSL doesn't seen to compute these in a
   * consistent way.
   */
  {
  static char* VALID_STRANDS[] = {
      "+", "-", "++", "+-", "-+", "--", NULL
  };
  int i, errCount = 0;
  int tBlockSizeMult = pslIsProtein(psl) ? 3 : 1;
  
  /* check strand value */
  for (i = 0; VALID_STRANDS[i] != NULL; i++)
      {
      if (strcmp(psl->strand, VALID_STRANDS[i]) == 0)
          break;
      }
  if (VALID_STRANDS[i] == NULL)
      chkError(pslDesc, out, psl, &errCount,
               "\tinvalid PSL strand: \"%s\"\n", psl->strand);
  
  /* check query */
  chkRanges(pslDesc, out, psl, psl->qName, "query", 'q', pslQStrand(psl), psl->qSize, psl->qStart, psl->qEnd,
            psl->qStarts, 1, &errCount);
  if ((opts & PSL_CHECK_IGNORE_INSERT_CNTS) == 0)
      chkInsertCounts(pslDesc, out, psl, psl->qName, 'q', psl->qStarts, psl->qNumInsert, psl->qBaseInsert, &errCount);
  
  /* check target */
  chkRanges(pslDesc, out, psl, psl->tName, "target", 't', pslTStrand(psl), psl->tSize, psl->tStart, psl->tEnd,
            psl->tStarts, tBlockSizeMult, &errCount);
  if ((opts & PSL_CHECK_IGNORE_INSERT_CNTS) == 0)
      chkInsertCounts(pslDesc, out, psl, psl->tName, 't', psl->tStarts, psl->tNumInsert, psl->tBaseInsert, &errCount);
  
  return errCount;
  }
  
  int pslCheck(char *pslDesc, FILE* out, struct psl* psl)
  /* Validate a PSL for consistency.  pslDesc is printed the error messages
   * to file out (open /dev/null to discard). Return count of errors. */
  {
  return pslCheck2(0, pslDesc, out, psl);
  }
  
  struct hash *readPslToBinKeeper(char *sizeFileName, char *pslFileName)
  /* read a list of psls and return results in hash of binKeeper structure for fast query*/
  {
  struct binKeeper *bk; 
  struct psl *psl;
  struct lineFile *sf = lineFileOpen(sizeFileName, TRUE);
  struct lineFile *pf = lineFileOpen(pslFileName , TRUE);
  struct hash *hash = newHash(0);
  char *chromRow[2];
  char *row[21] ;
  
  while (lineFileRow(sf, chromRow))
      {
      char *name = chromRow[0];
      int size = lineFileNeedNum(sf, chromRow, 1);
  
      if (hashLookup(hash, name) != NULL)
          warn("Duplicate %s, ignoring all but first\n", name);
      else
          {
          bk = binKeeperNew(0, size);
          assert(size > 1);
  	hashAdd(hash, name, bk);
          }
      }
  while (lineFileNextRow(pf, row, ArraySize(row)))
      {
      psl = pslLoad(row);
      bk = hashMustFindVal(hash, psl->tName);
      binKeeperAdd(bk, psl->tStart, psl->tEnd, psl);
      }
  lineFileClose(&pf);
  lineFileClose(&sf);
  return hash;
  }
  
  static boolean isDelChar(char c)
  /* is this a indel character? */
  {
  return (c == '-') || (c == '.') || (c == '=') || (c == '_');
  }
  
  static void trimAlignment(struct psl* psl, char** qStrPtr, char** tStrPtr,
                            int* aliSizePtr)
  /* remove leading or trailing indels from alignment */
  {
  char* qStr = *qStrPtr;
  char* tStr = *tStrPtr;
  int aliSize = *aliSizePtr;
  
  // skip lending indels
  while ((aliSize > 0) && (isDelChar(*qStr) || isDelChar(*tStr)))
      {
      if (!isDelChar(*qStr))
          psl->qStart++;
      else if (!isDelChar(*tStr))
          psl->tStart++;
      qStr++;
      tStr++;
      aliSize--;
      }
  
  // skip trailing indels
  while ((aliSize > 0) && (isDelChar(qStr[aliSize-1]) || isDelChar(tStr[aliSize-1])))
      {
      if (!isDelChar(qStr[aliSize-1]))
          psl->qEnd--;
      else if (!isDelChar(tStr[aliSize-1]))
          psl->tEnd--;
      aliSize--;
      }
  *qStrPtr = qStr;
  *tStrPtr = tStr;
  *aliSizePtr = aliSize;
  }
  
  static void addBlock(struct psl* psl, int qs, int qe, int ts, int te,
                       int *blockSpace)
  /* add a block to the psl, growing if necessary */
  {
  assert((qe-qs) == (te-ts));  /* same lengths? */
  if (psl->blockCount == *blockSpace)
      pslGrow(psl, blockSpace);
  psl->blockSizes[psl->blockCount] = qe - qs;
  psl->qStarts[psl->blockCount] = qs;
  psl->tStarts[psl->blockCount] = ts;
  psl->blockCount++;
  }
  
  static void accumCounts(struct psl *psl, char prevQ, char prevT,
                          char q, char t, unsigned options)
  /* accumulate block and base counts  */
  {
  if (!isDelChar(q) && !isDelChar(t))
      {
      /* aligned column. */
      char qu = toupper(q);
      char tu = toupper(t);
      if ((q == 'N') || (t == 'N'))
          psl->nCount++;
      else if (qu == tu)
          {
          if ((options & PSL_IS_SOFTMASK) && ((qu != q) || (tu != t)))
              psl->repMatch++;
          else
              psl->match++;
          }
      else
          psl->misMatch++;
      }
  else if (isDelChar(q) && !isDelChar(t))
      {
      /* target insert */
      psl->tBaseInsert++;
      if (!isDelChar(prevQ))
          psl->tNumInsert++;
      }
  else if (isDelChar(t) && !isDelChar(q))
      {
      /* query insert */
      psl->qBaseInsert++;
      if (!isDelChar(prevT))
          psl->qNumInsert++;
      }
  }
  
  struct psl* pslFromAlign(char *qName, int qSize, int qStart, int qEnd, char *qString,
                           char *tName, int tSize, int tStart, int tEnd, char *tString,
                           char* strand, unsigned options)
  /* Create a PSL from an alignment.  Options PSL_IS_SOFTMASK if lower case
   * bases indicate repeat masking.  Returns NULL if alignment is empty after
   * triming leading and trailing indels.*/
  {
  /* Note, the unit tests for these programs exercise this function:
   *   hg/embossToPsl
   *   hg/mouseStuff/axtToPsl
   *   hg/spideyToPsl
   *   utils/est2genomeToPsl
   */
  
  int blockSpace = 16;
  struct psl* psl = NULL;
  int aliSize = strlen(qString);
  boolean eitherInsert = FALSE;	/* True if either in insert state. */
  int i, qs,qe,ts,te;
  char prevQ = '\0',  prevT = '\0';
  AllocVar(psl);
  
  if (strlen(tString) != aliSize)
      errAbort("query and target alignment strings are different lengths");
  
  psl = pslNew(qName, qSize, qStart, qEnd, tName, tSize, tStart, tEnd,
               strand, blockSpace, 0);
  trimAlignment(psl, &qString, &tString, &aliSize);
  
  /* Don't create if either query or target is zero length after trim */
   if ((psl->qStart == psl->qEnd) || (psl->tStart == psl->tEnd))
       {
       pslFree(&psl);
       return NULL;
       }
  
  /* Get range of alignment in strand-specified coordinates */
  qs = psl->qStart;
  qe = psl->qEnd;
  if (strand[0] == '-')
      reverseIntRange(&qs, &qe, psl->qSize);
  ts = psl->tStart;
  te = psl->tEnd;
  if (strand[1] == '-')
      reverseIntRange(&ts, &te, psl->tSize);
  
  eitherInsert = FALSE;
  qe = qs;  /* current block coords */
  te = ts;
  for (i=0; i<aliSize; ++i)
      {
      char q = qString[i];
      char t = tString[i];
      if (isDelChar(q) && isDelChar(t))
          {
          continue; /* nothing in this column, just ignore it */
          }
      else if (isDelChar(q) || isDelChar(t))
          {
          /* insert in one of the columns */
  	if (!eitherInsert)
  	    {
              /* end of a block */
              addBlock(psl, qs, qe, ts, te, &blockSpace);
  	    eitherInsert = TRUE;
  	    }
  	if (!isDelChar(q))
              qe += 1;
  	if (!isDelChar(t))
              te += 1;
  	}
      else
          {
          /* columns aligned */
  	if (eitherInsert)
  	    {
              /* start new block */
  	    qs = qe;
  	    ts = te;
  	    eitherInsert = FALSE;
  	    }
  	qe += 1;
  	te += 1;
  	}
      accumCounts(psl, prevQ, prevT, q, t, options);
      prevQ = q; /* will not include skipped empty columns */
      prevT = t;
      }
  addBlock(psl, qs, qe, ts, te, &blockSpace);
  return psl;
  }
  
  struct psl* pslNew(char *qName, unsigned qSize, int qStart, int qEnd,
                     char *tName, unsigned tSize, int tStart, int tEnd,
                     char *strand, unsigned blockSpace, unsigned opts)
  /* create a new psl with space for the specified number of blocks allocated.
   * pslGrow maybe used to expand this space if needed.  Valid options are
   * PSL_XA_FORMAT. */
  {
  struct psl *psl;
  AllocVar(psl);
  assert(blockSpace > 0);
  psl->qName = cloneString(qName);
  psl->qSize = qSize;
  psl->qStart = qStart;
  psl->qEnd = qEnd;
  psl->tName = cloneString(tName);
  psl->tSize = tSize;
  psl->tStart = tStart;
  psl->tEnd = tEnd;
  strncpy(psl->strand, strand, 2);
  AllocArray(psl->blockSizes, blockSpace);
  AllocArray(psl->qStarts, blockSpace);
  AllocArray(psl->tStarts, blockSpace);
  if (opts & PSL_XA_FORMAT)
      {
      AllocArray(psl->qSequence, blockSpace);
      AllocArray(psl->tSequence, blockSpace);
      }
  return psl;
  }
  
  void pslGrow(struct psl *psl, int *blockSpacePtr)
  /* Increase memory allocated to a psl to hold more blocks.  blockSpacePtr
   * should point the the current maximum number of blocks and will be
   * updated to with the new amount of space. */
  {
  int blockSpace = *blockSpacePtr;
  int newSpace = 2 * blockSpace;
  ExpandArray(psl->blockSizes, blockSpace, newSpace);
  ExpandArray(psl->qStarts, blockSpace, newSpace);
  ExpandArray(psl->tStarts, blockSpace, newSpace);
  if (psl->qSequence != NULL)
      {
      ExpandArray(psl->qSequence, blockSpace, newSpace);
      ExpandArray(psl->tSequence, blockSpace, newSpace);
      }
  *blockSpacePtr = newSpace;
  }
  
  void pslComputeInsertCounts(struct psl *psl)
  /* compute numInsert and baseInsert fields from the blocks */
  {
  psl->qNumInsert = psl->qBaseInsert = 0;
  psl->tNumInsert = psl->tBaseInsert = 0;
  int iBlk;
  
  for (iBlk = 1; iBlk < psl->blockCount; iBlk++)
      {
      unsigned qGapSize = psl->qStarts[iBlk] - (psl->qStarts[iBlk-1]+psl->blockSizes[iBlk-1]);
      if (qGapSize != 0)
          {
          psl->qNumInsert++;
          psl->qBaseInsert += qGapSize;
          }
      unsigned tGapSize = psl->tStarts[iBlk] - (psl->tStarts[iBlk-1]+psl->blockSizes[iBlk-1]);
      if (tGapSize != 0)
          {
          psl->tNumInsert++;
          psl->tBaseInsert += tGapSize;
          }
      }
  }
  
  static boolean getNextCigarOp(char *startPtr, boolean reverse, char **ptr, char *op, int *size)
  /* gets one cigar op out of the CIGAR string.  Reverse the order if asked */
  {
  char *str = *ptr;
  
  if (str == NULL)
      return FALSE;
  
  if ((!reverse && (*str == 0)) || (reverse && (str == startPtr)))
      return FALSE;
  
  // between each cigar op there could be nothing, or a space, or a plus
  if (reverse)
      {
      char *end = str - 1;
      for(;*end ; end--)
  	{
  	if (isalpha(*end))
  	    break;
  	}
      str = end;
      *ptr = end;
      }
  else
      {
      char *end = str + 1;
      for(;*end ; end++)
  	{
  	if (! (isdigit(*end)  || (*end == ' ') || (*end == '+')))
  	    break;
  	}
  
      *ptr = end;
      }
  
  *op = *str++;
  *size = atoi(str);
  return TRUE;
  }
  
  struct psl* pslFromGff3Cigar(char *qName, int qSize, int qStart, int qEnd,
                               char *tName, int tSize, int tStart, int tEnd,
                               char* strand, char *cigar)
  /* create a PSL from a GFF3-style cigar formatted alignment */
  {
  // this is currently tested by the gff3ToPsl
  int blocksAlloced = 4;
  struct psl *psl = pslNew(qName, qSize, qStart, qEnd, tName, tSize, tStart, tEnd, strand, blocksAlloced, 0);
  
  char op;
  int size;
  int totalSize = 0;
  
  int qNext = qStart, qBlkEnd = qEnd;
  if (strand[0] == '-')
      reverseIntRange(&qNext, &qBlkEnd, qSize);
  int tNext = tStart, tBlkEnd = tEnd;
  if (strand[1] == '-')
      reverseIntRange(&tNext, &tBlkEnd, tSize);
  
  if (cigar == NULL)
      {
      // no cigar means one block
      size = qEnd - qStart;
      totalSize += size;
      psl->blockSizes[psl->blockCount] = size;
      psl->qStarts[psl->blockCount] = qNext;
      psl->tStarts[psl->blockCount] = tNext;
      psl->blockCount++;
      tNext += size;
      qNext += size;
      }
  else
      {
      if (strand[0] == '-' && strand[1] == '-')
          errAbort("GFF3 spec is vague about Gap when both strands are '-'; not implemented yet.");
      char cigarSpec[strlen(cigar)+1];  // copy since parsing is destructive
      safecpy(cigarSpec, sizeof cigarSpec, cigar);
      char *cigarNext = cigarSpec;
      while(getNextCigarOp(cigarSpec, FALSE, &cigarNext, &op, &size))
  	{
  	switch (op)
  	    {
  	    case 'M': // match or mismatch (gapless aligned block)
  		if (psl->blockCount == blocksAlloced)
  		    pslGrow(psl, &blocksAlloced);
  
  		totalSize += size;
  		psl->blockSizes[psl->blockCount] = size;
  		psl->qStarts[psl->blockCount] = qNext;
  		psl->tStarts[psl->blockCount] = tNext;
  		psl->blockCount++;
  		tNext += size;
  		qNext += size;
  		break;
  	    case 'I': // inserted in target
  		tNext += size;
  		break;
  	    case 'D': // deleted from target
  		qNext += size;
  		break;
  	    
  	    default:
  		errAbort("unrecognized CIGAR op %c in %s", op, cigar);
  	    }
  	}
      }
  
  /* CIGARs starting/ending with indels require adjusting of query/target ranges,
   * as PSL starts/ends with matches */
  psl->qStart = psl->qStarts[0];
  psl->qEnd = pslQEnd(psl, psl->blockCount-1);
  if (strand[0] == '-')
      reverseIntRange(&psl->qStart, &psl->qEnd, qSize);
  psl->tStart = psl->tStarts[0];
  psl->tEnd = pslTEnd(psl, psl->blockCount-1);
  if (strand[1] == '-')
      reverseIntRange(&psl->tStart, &psl->tEnd, tSize);
  
  /* sanity check */
  if (qNext != qBlkEnd)
      errAbort("CIGAR query length does not match specified query range %s:%d-%d", qName, qStart, qEnd);
  if (tNext != tBlkEnd)
      errAbort("CIGAR target length does not match specified target range %s:%d-%d", tName, tStart, tEnd);
  psl->match = totalSize;
  pslComputeInsertCounts(psl);
  return psl;
  }
  
  int pslRangeTreeOverlap(struct psl *psl, struct rbTree *rangeTree)
  /* Return amount that psl overlaps (on target side) with rangeTree. */
  {
  int i;
  int overlap = 0;
  boolean isRc = (psl->strand[1] == '-');
  for (i=0; i<psl->blockCount; ++i)
      {
      int start = psl->tStarts[i];
      int end = start + psl->blockSizes[i];
      if (isRc)
          reverseIntRange(&start, &end, psl->tSize);
      overlap += rangeTreeOverlapSize(rangeTree, start, end);
      }
  return overlap;
  }
  
  float pslIdent(struct psl *psl)
  /* computer fraction identity */
  {
  float aligned = psl->match + psl->misMatch + psl->repMatch;
  if (aligned == 0.0)
      return 0.0;
  else
      return ((float)(psl->match + psl->repMatch))/aligned;
  }
  
  float pslQueryAligned(struct psl *psl)
  /* compute fraction of query that was aligned */
  {
  float aligned = psl->match + psl->misMatch + psl->repMatch;
  return aligned/(float)psl->qSize;
  }
  
  struct psl* pslClone(struct psl *psl)
  /* clone a psl */
  {
  struct psl* pslCp = pslNew(psl->qName, psl->qSize, psl->qStart, psl->qEnd,
                             psl->tName, psl->tSize, psl->tStart, psl->tEnd,
                             psl->strand, psl->blockCount,
                             ((psl->tSequence != NULL) ? PSL_XA_FORMAT : 0));
  pslCp->match = psl->match;
  pslCp->misMatch = psl->misMatch;
  pslCp->repMatch = psl->repMatch;
  pslCp->nCount = psl->nCount;
  pslCp->qNumInsert = psl->qNumInsert;
  pslCp->qBaseInsert = psl->qBaseInsert;
  pslCp->tNumInsert = psl->tNumInsert;
  pslCp->tBaseInsert = psl->tBaseInsert;
  int iBlk;
  for (iBlk = 0; iBlk < psl->blockCount; iBlk++)
      {
      pslCp->blockSizes[iBlk] = psl->blockSizes[iBlk];
      pslCp->qStarts[iBlk] = psl->qStarts[iBlk];
      pslCp->tStarts[iBlk] = psl->tStarts[iBlk];
      if (psl->qSequence != NULL)
          pslCp->qSequence[iBlk] = cloneString(psl->qSequence[iBlk]);
      if (psl->tSequence != NULL)
          pslCp->tSequence[iBlk] = cloneString(psl->tSequence[iBlk]);
      pslCp->blockCount++;
      }
  return pslCp;
  }
  
  int cmpChrom(char *a, char *b)
  /* Compare two chromosomes. */
  {
  return cmpStringsWithEmbeddedNumbers(a, b);
  }
  
  
  int pslCmpTargetScore(const void *va, const void *vb)
  /* Compare to sort based on target then score. */
  {
  const struct psl *a = *((struct psl **)va);
  const struct psl *b = *((struct psl **)vb);
  int diff = cmpChrom(a->tName, b->tName);
  if (diff == 0)
      diff = pslScore(b) - pslScore(a);
  return diff;
  }
  
  int pslCmpTargetStart(const void *va, const void *vb)
  /* Compare to sort based on target start. */
  {
  const struct psl *a = *((struct psl **)va);
  const struct psl *b = *((struct psl **)vb);
  int diff = cmpChrom(a->tName, b->tName);
  if (diff == 0)
      diff = a->tStart - b->tStart;
  return diff;
  }
  
  char *pslSortList[] = {"query,score", "query,start", "chrom,score", "chrom,start", "score"};
  
  void pslSortListByVar(struct psl **pslList, char *sort)
  /* Sort a list of psls using the method definied in the sort string. */
  {
  if (sameString(sort, "query,start"))
      {
      slSort(pslList, pslCmpQuery);
      }
  else if (sameString(sort, "query,score"))
      {
      slSort(pslList, pslCmpQueryScore);
      }
  else if (sameString(sort, "score"))
      {
      slSort(pslList, pslCmpScore);
      }
  else if (sameString(sort, "chrom,start"))
      {
      slSort(pslList, pslCmpTargetStart);
      }
  else if (sameString(sort, "chrom,score"))
      {
      slSort(pslList, pslCmpTargetScore);
      }
  else
      {
      slSort(pslList, pslCmpQueryScore);
      }
  }
  
  void pslRemoveFrameShifts(struct psl *psl)
  /* Remove any frameshits if present. Changes in place, doesn't update statistics in first nine fields. */
  {
  unsigned prevEnd = psl->tStarts[0];
  int ii;
  
  for(ii = 0; ii < psl->blockCount; ii++)
      {
      int diff = prevEnd - psl->tStarts[ii];
      if (diff > 0)
          {
          if (diff > psl->blockSizes[ii])
              errAbort("frame shift (%d) larger than block size (%d)", diff, psl->blockSizes[ii]);
          psl->blockSizes[ii] -= diff;
          psl->tStarts[ii] += diff;
          psl->qStarts[ii] += diff;
          }
      prevEnd = psl->tStarts[ii] + psl->blockSizes[ii];
      }
  }
  