f2b126d2fe65e705b4f6554f1c0c3b8dfb6d0a83
angie
  Wed May 31 12:37:39 2023 -0700
Check our reconstructed ancestral base value vs. the tree's par (parental/parsimonious?) base value.

diff --git src/hg/hgPhyloPlace/treeToAuspiceJson.c src/hg/hgPhyloPlace/treeToAuspiceJson.c
index 29efa08..b43eb40 100644
--- src/hg/hgPhyloPlace/treeToAuspiceJson.c
+++ src/hg/hgPhyloPlace/treeToAuspiceJson.c
@@ -390,50 +390,60 @@
                                           struct singleNucChange *ancestorMuts,
                                           struct geneInfo *gi)
 /* Given a list of vpTx from the same codon, combine their changes with inherited mutations
  * in the same codon to get the amino acid change at this node.
  * Note: this assumes there is no UTR in transcript, only CDS.  True so far for pathogens... */
 {
 struct slName *aaChange = NULL;
 if (codonVpTxList != NULL)
     {
     struct vpTx *vpTx = codonVpTxList;
     int firstAaStart = vpTx->start.txOffset / 3;
     int codonStart = firstAaStart * 3;
     int codonOffset = vpTx->start.txOffset - codonStart;
     char oldCodon[4];
     safencpy(oldCodon, sizeof oldCodon, gi->txSeq->dna + codonStart, 3);
+    touppers(oldCodon);
     int codonStartG = vpTx->start.gOffset - codonOffset;
     struct singleNucChange *anc;
     for (anc = ancestorMuts;  anc != NULL;  anc = anc->next)
         {
         int ancCodonOffset = anc->chromStart - codonStartG;
         if (ancCodonOffset >= 0 && ancCodonOffset < 3)
+            {
+            if (anc->parBase != oldCodon[ancCodonOffset])
+                errAbort("codonVpTxToAaChange: expected parBase for ancestral mutation at %d "
+                         "(%s codon %d offset %d) to be '%c' but it's '%c'",
+                         anc->chromStart+1, gi->psl->qName, firstAaStart, ancCodonOffset,
+                         oldCodon[ancCodonOffset], anc->parBase);
             oldCodon[ancCodonOffset] = maybeComplement(anc->newBase, gi->psl);
             }
+        }
     char oldAa = lookupCodon(oldCodon);
     if (oldAa == '\0')
         oldAa = '*';
     char newCodon[4];
     safecpy(newCodon, sizeof newCodon, oldCodon);
     for (vpTx = codonVpTxList;  vpTx != NULL;  vpTx = vpTx->next)
         {
         int aaStart = vpTx->start.txOffset / 3;
         if (aaStart != firstAaStart)
             errAbort("codonVpTxToAaChange: program error: firstAaStart %d != aaStart %d",
                      firstAaStart, aaStart);
         int codonOffset = vpTx->start.txOffset - codonStart;
+        // vpTx->txRef[0] is always the reference base, not like singleNucChange parBase,
+        // so we can't compare it to expected value as we could for ancMuts above.
         newCodon[codonOffset] = vpTx->txAlt[0];
         }
     char newAa = lookupCodon(newCodon);
     if (newAa == '\0')
         newAa = '*';
     if (newAa != oldAa)
         {
         char aaChangeStr[32];
         safef(aaChangeStr, sizeof aaChangeStr, "%c%d%c", oldAa, firstAaStart+1, newAa);
         aaChange = slNameNew(aaChangeStr);
         }
     }
 return aaChange;
 }