b269cf8292be3c573bb92058fbd6a3df5cf385c9
hiram
  Thu Sep 24 15:44:03 2026 -0700
translation of Chenxi Zhou code from Gene Myers FASTGA PAFtoPSL.c code to kent C code refs #34360

diff --git src/utils/pafToPsl/pafToPsl.c src/utils/pafToPsl/pafToPsl.c
new file mode 100644
index 00000000000..bceab7ac11f
--- /dev/null
+++ src/utils/pafToPsl/pafToPsl.c
@@ -0,0 +1,327 @@
+/* pafToPsl - converting paf output from minimap2 into UCSC psl format.
+ * This is a conversion on permission from author: Chenxi Zhou from code
+ * in the Gene Myers project: fastGA/FASTGA utility PAFtoPSL
+ *   https://github.com/thegenemyers/FASTGA
+ * to function in the kent C code environment.
+ */
+
+#include "common.h"
+#include "linefile.h"
+#include "options.h"
+#include "psl.h"
+#include "sqlNum.h"
+#include "hash.h"
+#include "obscure.h"
+#include <ctype.h>
+
+void usage()
+/* Explain usage and exit. */
+{
+errAbort(
+  "pafToPsl - converting paf output from minimap2 into UCSC psl format\n"
+  "usage:\n"
+  "   pafToPsl in.paf out.psl\n"
+  "options:\n"
+  "   -cigarTag=tag  SAM-style tag that carries the alignment CIGAR string\n"
+  "                  (default cg:Z:).  The CIGAR is made of M/=/X, I and D\n"
+  "                  operators, as minimap2 writes it.\n"
+  "   -tSizes=file.chrom.sizes  Tab-separated <name><size> file.  When given,\n"
+  "                  each record's target size is taken from here instead of\n"
+  "                  the PAF's own tSize column; a target name missing from\n"
+  "                  this file causes that record to be skipped with a\n"
+  "                  warning, and a size that disagrees with the PAF is\n"
+  "                  also warned about (this file's value wins).\n"
+  "   -qSizes=file.chrom.sizes  Same as -tSizes, for the query side.\n"
+  "out.psl may be \"stdout\" to write to standard output.\n"
+  "### Thank you to Chenxi Zhou for permission to translate the original  ###\n"
+  "### PAFtoPSL code from the Gene Myers project: FASTGA                  ###\n"
+  "###     https://github.com/thegenemyers/FASTGA                         ###\n"
+  "### into this kent C code implementation.                              ###"
+  );
+}
+
+static char *cigarTag = "cg:Z:";
+static char *tSizesFile = NULL;
+static char *qSizesFile = NULL;
+
+/* Command line validation table. */
+static struct optionSpec options[] = {
+   {"cigarTag", OPTION_STRING},
+   {"tSizes", OPTION_STRING},
+   {"qSizes", OPTION_STRING},
+   {NULL, 0},
+};
+
+static void pushBlock(struct psl *psl, int *blockSpace, int qPos, int tPos, int size)
+/* Append one ungapped alignment block to psl, growing the block arrays if
+ * needed. */
+{
+if (psl->blockCount >= *blockSpace)
+    pslGrow(psl, blockSpace);
+psl->qStarts[psl->blockCount] = qPos;
+psl->tStarts[psl->blockCount] = tPos;
+psl->blockSizes[psl->blockCount] = size;
+psl->blockCount++;
+}
+
+static boolean cigarToPsl(struct psl *psl, int *blockSpace, char *cigar,
+	char *fileName, int lineIx)
+/* Parse a CIGAR string (as found in a PAF cg:Z: tag) into psl's alignment
+ * blocks, and fill in the fields that are derived from it: the trimmed
+ * qStart/qEnd/tStart/tEnd, the insert counts, misMatch and nCount.
+ * psl->match must already be set from the PAF record's own match count.
+ * Returns FALSE, after printing a warning, if the CIGAR is inconsistent
+ * with psl's alignment range. */
+{
+int qNumInsert = 0, qBaseInsert = 0, tNumInsert = 0, tBaseInsert = 0;
+int qPos = 0, tPos = 0, lens = 0, insl = 0, insr = 0, opLen = 0;
+char *c = cigar;
+char op, prevOp = '\0';
+int i;
+
+psl->blockCount = 0;
+while (*c != '\0')
+    {
+    opLen = 0;
+    while (isdigit((unsigned char)*c))
+        opLen = 10*opLen + (*c++ - '0');
+    if (opLen == 0)
+	{
+	warn("%s:%d: CIGAR operator length is zero", fileName, lineIx);
+	return FALSE;
+	}
+    op = *c++;
+    switch (op)
+        {
+	case 'M':
+	case 'X':
+	case '=':
+	    qPos += opLen;
+	    tPos += opLen;
+	    lens += opLen;
+	    break;
+	case 'I':
+	    if (prevOp == '\0')	// leading insertion
+		insl = opLen;
+	    else
+		{
+		pushBlock(psl, blockSpace, qPos-lens, tPos-lens, lens);
+		lens = 0;
+		}
+	    qNumInsert += 1;
+	    qBaseInsert += opLen;
+	    qPos += opLen;
+	    break;
+	case 'D':
+	    if (prevOp == '\0')	// leading deletion
+		insl = -opLen;
+	    else
+		{
+		pushBlock(psl, blockSpace, qPos-lens, tPos-lens, lens);
+		lens = 0;
+		}
+	    tNumInsert += 1;
+	    tBaseInsert += opLen;
+	    tPos += opLen;
+	    break;
+	default:
+	    warn("%s:%d: invalid CIGAR operator '%c'", fileName, lineIx, op);
+	    return FALSE;
+	}
+    prevOp = op;
+    }
+
+// the last CIGAR operator
+if (prevOp == 'I')		// trailing insertion
+    insr = opLen;
+else if (prevOp == 'D')	// trailing deletion
+    insr = -opLen;
+else
+    pushBlock(psl, blockSpace, qPos-lens, tPos-lens, lens);
+
+if (qPos != psl->qEnd - psl->qStart)
+    {
+    warn("%s:%d: CIGAR length does not match alignment length (query): %d != %d",
+	fileName, lineIx, psl->qEnd - psl->qStart, qPos);
+    return FALSE;
+    }
+if (tPos != psl->tEnd - psl->tStart)
+    {
+    warn("%s:%d: CIGAR length does not match alignment length (target): %d != %d",
+	fileName, lineIx, psl->tEnd - psl->tStart, tPos);
+    return FALSE;
+    }
+
+// handle leading and trailing indels: PSL blocks must start and end on a match
+if (insl > 0)
+    { qNumInsert -= 1; qBaseInsert -= insl; psl->qStart += insl; }
+else if (insl < 0)
+    { tNumInsert -= 1; tBaseInsert += insl; psl->tStart -= insl; }
+
+if (insr > 0)
+    { qNumInsert -= 1; qBaseInsert -= insr; psl->qEnd -= insr; }
+else if (insr < 0)
+    { tNumInsert -= 1; tBaseInsert += insr; psl->tEnd += insr; }
+
+// shift target block positions to absolute target coordinates
+for (i = 0; i < psl->blockCount; i++)
+    psl->tStarts[i] += psl->tStart;
+
+// shift query block positions to absolute coordinates, taking the query
+// strand into account
+if (psl->strand[0] == '-')
+    for (i = 0; i < psl->blockCount; i++)
+	psl->qStarts[i] = psl->qSize - psl->qEnd + psl->qStarts[i];
+else
+    for (i = 0; i < psl->blockCount; i++)
+	psl->qStarts[i] += psl->qStart;
+
+psl->qNumInsert = qNumInsert;
+psl->qBaseInsert = qBaseInsert;
+psl->tNumInsert = tNumInsert;
+psl->tBaseInsert = tBaseInsert;
+
+int misMatch = psl->qEnd - psl->qStart - qBaseInsert - (int)psl->match;
+if (misMatch < 0)
+    {
+    warn("%s:%d: negative misMatch value: %d", fileName, lineIx, misMatch);
+    return FALSE;
+    }
+psl->misMatch = misMatch;
+
+lens = 0;
+for (i = 0; i < psl->blockCount; i++)
+    lens += psl->blockSizes[i];
+int nCount = lens - (int)psl->match - (int)psl->misMatch - (int)psl->repMatch;
+if (nCount < 0)
+    {
+    warn("%s:%d: negative nCount value: %d", fileName, lineIx, nCount);
+    return FALSE;
+    }
+psl->nCount = nCount;
+
+return TRUE;
+}
+
+static void pafToPsl(char *inFile, char *outFile)
+/* pafToPsl - converting paf output from minimap2 into UCSC psl format. */
+{
+struct lineFile *lf = lineFileOpen(inFile, TRUE);
+FILE *f = mustOpen(outFile, "w");
+int tagLen = strlen(cigarTag);
+struct hash *tSizeHash = tSizesFile ? hashNameIntFile(tSizesFile) : NULL;
+struct hash *qSizeHash = qSizesFile ? hashNameIntFile(qSizesFile) : NULL;
+char *line;
+int lineSize;
+
+while (lineFileNext(lf, &line, &lineSize))
+    {
+    char *s = line;
+    char *fields[11];
+    int i;
+
+    for (i = 0; i < ArraySize(fields); i++)
+	{
+	fields[i] = nextWord(&s);
+	if (fields[i] == NULL)
+	    break;
+	}
+    if (i < ArraySize(fields))
+	{
+	warn("%s:%d: PAF line has fewer than %d fields, skipping",
+	    lf->fileName, lf->lineIx, (int)ArraySize(fields));
+	continue;
+	}
+
+    char *cigar = NULL;
+    char *tag;
+    while ((tag = nextWord(&s)) != NULL)
+	{
+	if ((int)strlen(tag) > tagLen && startsWith(cigarTag, tag))
+	    {
+	    cigar = tag + tagLen;
+	    break;
+	    }
+	}
+    if (cigar == NULL)
+	{
+	warn("%s:%d: PAF line is missing a %s CIGAR tag, skipping",
+	    lf->fileName, lf->lineIx, cigarTag);
+	continue;
+	}
+
+    char strand[2];
+    strand[0] = fields[4][0];
+    strand[1] = '\0';
+    if (strand[0] != '+' && strand[0] != '-')
+	{
+	warn("%s:%d: invalid PAF strand '%c', skipping",
+	    lf->fileName, lf->lineIx, strand[0]);
+	continue;
+	}
+
+    unsigned qSize = sqlUnsigned(fields[1]);
+    if (qSizeHash != NULL)
+	{
+	int sizesVal = hashIntValDefault(qSizeHash, fields[0], -1);
+	if (sizesVal < 0)
+	    {
+	    warn("%s:%d: query '%s' not found in %s, skipping",
+		lf->fileName, lf->lineIx, fields[0], qSizesFile);
+	    continue;
+	    }
+	if ((unsigned)sizesVal != qSize)
+	    warn("%s:%d: query '%s' size %u in PAF does not match %d in %s, using %s",
+		lf->fileName, lf->lineIx, fields[0], qSize, sizesVal, qSizesFile, qSizesFile);
+	qSize = (unsigned)sizesVal;
+	}
+
+    unsigned tSize = sqlUnsigned(fields[6]);
+    if (tSizeHash != NULL)
+	{
+	int sizesVal = hashIntValDefault(tSizeHash, fields[5], -1);
+	if (sizesVal < 0)
+	    {
+	    warn("%s:%d: target '%s' not found in %s, skipping",
+		lf->fileName, lf->lineIx, fields[5], tSizesFile);
+	    continue;
+	    }
+	if ((unsigned)sizesVal != tSize)
+	    warn("%s:%d: target '%s' size %u in PAF does not match %d in %s, using %s",
+		lf->fileName, lf->lineIx, fields[5], tSize, sizesVal, tSizesFile, tSizesFile);
+	tSize = (unsigned)sizesVal;
+	}
+
+    int blockSpace = 16;
+    struct psl *psl = pslNew(fields[0], qSize,
+	sqlSigned(fields[2]), sqlSigned(fields[3]),
+	fields[5], tSize,
+	sqlSigned(fields[7]), sqlSigned(fields[8]),
+	strand, blockSpace, 0);
+    psl->match = sqlUnsigned(fields[9]);
+    // fields[10], the PAF "number of minimizers"/alignment block length
+    // field, is not used -- as in the original converter.
+
+    if (cigarToPsl(psl, &blockSpace, cigar, lf->fileName, lf->lineIx))
+	pslTabOut(psl, f);
+    pslFree(&psl);
+    }
+
+carefulClose(&f);
+lineFileClose(&lf);
+}
+
+int main(int argc, char *argv[])
+/* Process command line. */
+{
+optionInit(&argc, argv, options);
+if (argc != 3)
+    usage();
+cigarTag = optionVal("cigarTag", cigarTag);
+if (strlen(cigarTag) != 5 || cigarTag[2] != ':' || cigarTag[4] != ':')
+    errAbort("-cigarTag must look like xx:Z: (got '%s')", cigarTag);
+tSizesFile = optionVal("tSizes", NULL);
+qSizesFile = optionVal("qSizes", NULL);
+pafToPsl(argv[1], argv[2]);
+return 0;
+}