ad81d8de49a9c061fb50958f72a017ae1436497d angie Fri Aug 28 08:17:19 2026 -0700 Fix regex for rerooting pango.clade-mutations.tsv to also match B's empty mutation list. Bug caught by Claude automated review of 0b09514. Fortunately masked in latest pangolin-data release because the lineage B annotation was preserved when rerooting the tree and did not need re-annotating. diff --git src/hg/utils/otto/sarscov2phylo/updateLineageTreePb.sh src/hg/utils/otto/sarscov2phylo/updateLineageTreePb.sh index 01396ed5abe..f58c25a9481 100755 --- src/hg/utils/otto/sarscov2phylo/updateLineageTreePb.sh +++ src/hg/utils/otto/sarscov2phylo/updateLineageTreePb.sh @@ -1,195 +1,195 @@ #!/bin/bash source ~/.bashrc set -beEu -x -o pipefail # Do not modify this script, modify the source tree copy: # kent/src/hg/utils/otto/sarscov2phylo/updateLineageTreePb.sh usage() { echo "usage: $0 buildDate [tree.pb.gz]" } if (( $# != 1 && $# != 2 )); then usage exit 1 fi buildDate=$1 if [ $# == 2 ]; then startingTree=$2 else startingTree=gisaidAndPublic.$buildDate.masked.pb.gz fi ottoDir=/hive/data/outside/otto/sarscov2phylo usherDir=~angie/github/usher matUtils=$usherDir/build/matUtils today=$(date +%F) cd $ottoDir/$buildDate # Remove sequences that have two or more reversions relative to their assigned clade/lineage. $matUtils summary -i $startingTree --node-stats node-stats # Until BA.2.3.22 gets more samples that don't have bogus reversion on 25000 and 26577, # exempt some samples; also, most of JP.1 has rev on 27383,27384 and 2 of LF.7.6.5's 4 # samples have revs. cat > pruneRevsExemptions < 1 {print $1;}' node-stats \ | grep -vFwf pruneRevsExemptions \ > pruneRevs $matUtils extract -i $startingTree \ -p -s pruneRevs -O -o gisaidAndPublic.$buildDate.masked.pruneRevs.pb.gz # Get node ID for root of lineage A, used as reference/root by Pangolin: $matUtils extract -i gisaidAndPublic.$buildDate.masked.pruneRevs.pb.gz -C clade-paths.prunedRevs lineageARoot=$(grep ^A$'\t' clade-paths.prunedRevs | cut -f 2) # Reroot protobuf to lineage A and restrict to low mutation density (highly supported nodes): $matUtils extract -i gisaidAndPublic.$buildDate.masked.pruneRevs.pb.gz \ --reroot $lineageARoot \ --max-mutation-density 2 \ -O -o gisaidAndPublic.$buildDate.masked.reroot.pb.gz # Reroot pango.clade-mutations.tsv -sed -re 's/\t([A-Z][0-9]+[A-Z])/\tT8782C > C28144T > \1/;' $scriptDir/pango.clade-mutations.tsv \ +sed -re 's/\t($|[A-Z][0-9]+[A-Z])/\tT8782C > C28144T > \1/;' $scriptDir/pango.clade-mutations.tsv \ > pango.clade-mutations.reroot.tsv # Mask additional bases at the beginning and end of the genome that pangolin masks after # aligning input sequences. set +x for ((i=56; $i <= 265; i++)); do echo -e "N${i}N" done > maskPangoEnds for ((i=29674; $i < 29804; i++)); do echo -e "N${i}N" done >> maskPangoEnds set -x $matUtils mask -i gisaidAndPublic.$buildDate.masked.reroot.pb.gz -c \ -m maskPangoEnds -o gisaidAndPublic.$buildDate.masked.reroot.pangoMasked.pb.gz # Preserve lineage annotations that survived rerooting and pango-masking $matUtils extract -i gisaidAndPublic.$buildDate.masked.reroot.pangoMasked.pb.gz \ -C clade-paths.reroot.pangoMasked tail -n+2 clade-paths.reroot.pangoMasked \ | grep '^[A-Za-z]' \ | cut -f 1,3 > lineageToPath.reroot.pangoMasked # Assign updated lineages on the rerooted & pango-masked tree, pango-only for pangolin: time $matUtils annotate -T 50 \ -i gisaidAndPublic.$buildDate.masked.reroot.pangoMasked.pb.gz \ -P lineageToPath.reroot.pangoMasked \ -M pango.clade-mutations.reroot.tsv \ -l \ -c lineageToName \ -f 0.95 \ -u mutations.pangoOnly \ -D details.pangoOnly \ -o gisaidAndPublic.$buildDate.masked.reroot.pangoOnly.pb.gz \ >& annotate.pangoOnly.out set +o pipefail grep 'Could not' annotate.pangoOnly.out | cat grep skip annotate.pangoOnly.out | cat set -o pipefail # Make candidate subtrees. Remove BA.5.1_no29666 because it's always a problematic mini-BA.5.1; # then reassign all extra labels (for potential recombinants and sublineages) to valid Pango # labels so we don't get a regex failure when running pangolin. mkdir -p /hive/users/angie/lineageTreeUpdate.$today cd /hive/users/angie/lineageTreeUpdate.$today for i in 0 1 2 3 4 5 6 7 8 9; do echo test.50.$i $matUtils extract -i $ottoDir/$buildDate/gisaidAndPublic.$buildDate.masked.reroot.pangoOnly.pb.gz \ -r 50 -o test.50.$i.pb $matUtils mask -i test.50.$i.pb -S -o test.50.$i.simp.pb $matUtils extract -i test.50.$i.simp.pb -p -c BA.5.1_no29666 -o test.50.$i.simp.tmp.pb $matUtils extract -i test.50.$i.simp.tmp.pb -C clade-paths.$i.tmp tail -n+2 clade-paths.$i.tmp \ | cut -f 1,2 \ | sed -re 's/miscBA[A-Za-z0-9]+/BA.2/; s/miscDeltaBA1[A-Za-z0-9]+/BA.1/; s/miscBA.5.2CJ.1/CJ.1/; s/miscBA.5BA.2.75/BA.5/; s/proposed422/BA.1/; s/proposed437/BA.1/; s/proposed439/BA.1.1/; s/proposed441/BA.1/; s/proposed446/BA.1.1/; s/proposed455/BA.1/; s/proposed467/BA.2/; s/proposed482/BA.2/; s/proposed882/BA.2/; s/proposed885/BA.2/; s/proposed911/BA.5/; s/proposed1006/BA.2/; s/proposed1137/BA.5/; s/([A-Z\.0-9]+)_[A-Za-z0-9_]+/\1/;' \ > cladeNodes.reassign $matUtils annotate -i test.50.$i.simp.tmp.pb -l -C cladeNodes.reassign -o test.50.$i.simp.pb rm test.50.$i.simp.tmp.pb clade-paths.$i.tmp cladeNodes.reassign done # Run pangolin on each subtree # About 1 minute per job conda activate pangolin for i in 0 1 2 3 4 5 6 7 8 9; do echo test.50.$i.simp time pangolin -t 80 --usher-tree test.50.$i.simp.pb --skip-scorpio \ --skip-designation-cache --no-temp --outdir subset_10000_0.pusher.test.50.$i.simp.out \ ../pangolin_eval/subset_10000_0.fa done # Summarize results for i in 0 1 2 3 4 5 6 7 8 9; do echo test.50.$i.simp tail -n+2 subset_10000_0.pusher.test.50.$i.simp.out/lineage_report.csv \ | awk -F, '{print $1 "\t" $2;}' \ | sort \ > subset_10000_0.pusher.test.50.$i.simp join -t$'\t' ../pangolin_eval/subset_10000_0.cogNameToLin subset_10000_0.pusher.test.50.$i.simp \ | tawk '$2 != $3' \ > subset_10000_0.pusher.test.50.$i.simp.diff cut -f 2,3 subset_10000_0.pusher.test.50.$i.simp.diff \ | sort | uniq -c | sort -nr > subset_10000_0.pusher.test.50.$i.simp.diffrank done wc -l subset_10000_0.*.diff | sort -n # Pick the subtree with the fewest differences set +o pipefail bestTree=$(wc -l subset_10000_0.*.diff | sort -n | head -1 \ | sed -re 's/.*pusher\.test\.50\.//; s/\.diff//') set -o pipefail echo $bestTree # Which lineages have the most differences? head subset_10000_0.pusher.test.50.$bestTree.diffrank