4efa4aebadc39cb64a256c3489035b00bd2abdc6
angie
  Wed Sep 8 12:06:50 2021 -0700
Remove dependency on reference.fa from aligning new sequences; add tree-annotated Nextstrain clade and Pango lineage to metadata.

diff --git src/hg/utils/otto/sarscov2phylo/extractPublicTree.sh src/hg/utils/otto/sarscov2phylo/extractPublicTree.sh
index 5a8a750..b1fdda1 100755
--- src/hg/utils/otto/sarscov2phylo/extractPublicTree.sh
+++ src/hg/utils/otto/sarscov2phylo/extractPublicTree.sh
@@ -16,32 +16,36 @@
 
 today=$1
 ottoDir=/hive/data/outside/otto/sarscov2phylo
 cncbDir=$ottoDir/cncb.latest
 
 scriptDir=$(dirname "${BASH_SOURCE[0]}")
 source $scriptDir/util.sh
 
 usherDir=~angie/github/usher
 matUtils=$usherDir/build/matUtils
 
 cd $ottoDir/$today
 
 # Extract public samples from tree
 grep -v EPI_ISL_ samples.$today > newPublicNames
+# Dunno why, but when I tried using -s together with the filtering params, it ran for 3 hours
+# and I killed it -- stuck in a loop?  Run two commands:
 $matUtils extract -i gisaidAndPublic.$today.masked.pb \
     -s newPublicNames \
+    -o public-$today.all.masked.preTrim.pb
+$matUtils extract -i public-$today.all.masked.preTrim.pb \
     --max-parsimony 20 \
     --max-branch-length 30 \
     -O -o public-$today.all.masked.pb
 
 # Add nextclade annotations to public protobuf
 if [ -s cladeToName ]; then
     # Use combined tree's clade assignments to annotate clades on public tree
     grep -v EPI_ISL cladeToName > cladeToPublicName
     time $matUtils annotate -T 50 \
         -l \
         -i public-$today.all.masked.pb \
         -c cladeToPublicName \
         -f 0.95 \
         -D details.nextclade.public \
         -o public-$today.all.masked.nextclade.pb \
@@ -84,37 +88,38 @@
 rm public-$today.all.masked.pb
 ln -f public-$today.all.masked.nextclade.pangolin.pb public-$today.all.masked.pb
 
 cncbDate=$(ls -l $cncbDir | sed -re 's/.*cncb\.([0-9]{4}-[0-9][0-9]-[0-9][0-9]).*/\1/')
 echo "sarscov2phylo release 13-11-20; NCBI and COG-UK sequences downloaded $today; CNCB sequences downloaded $cncbDate" \
     > version.txt
 
 $matUtils extract -i public-$today.all.masked.pb -u samples.public.$today
 sampleCountComma=$(echo $(wc -l < samples.public.$today) \
                    | sed -re 's/([0-9]+)([0-9]{3})$/\1,\2/; s/([0-9]+)([0-9]{3},[0-9]{3})$/\1,\2/;')
 echo "$sampleCountComma genomes from GenBank, COG-UK and CNCB ($today); sarscov2phylo 13-11-20 tree with newer sequences added by UShER" \
     > hgPhyloPlace.description.txt
 
 # Make Taxodium-formatted protobuf for display
 zcat /hive/data/genomes/wuhCor1/goldenPath/bigZips/genes/ncbiGenes.gtf.gz > ncbiGenes.gtf
+zcat /hive/data/genomes/wuhCor1/wuhCor1.fa.gz > wuhCor1.fa
 zcat public-$today.metadata.tsv.gz > metadata.tmp.tsv
 time $matUtils extract -i public-$today.all.masked.pb \
-    -f reference.fa \
+    -f wuhCor1.fa \
     -g ncbiGenes.gtf \
     -M metadata.tmp.tsv \
     --write-taxodium public-$today.all.masked.taxodium.pb
-rm metadata.tmp.tsv
+rm metadata.tmp.tsv wuhCor1.fa
 gzip public-$today.all.masked.taxodium.pb
 
 # Link to public trees download directory hierarchy
 archiveRoot=/hive/users/angie/publicTrees
 read y m d < <(echo $today | sed -re 's/-/ /g')
 archive=$archiveRoot/$y/$m/$d
 mkdir -p $archive
 gzip -c public-$today.all.nwk > $archive/public-$today.all.nwk.gz
 ln -f `pwd`/public-$today.all.masked.{pb,vcf.gz} $archive/
 gzip -c public-$today.all.masked.pb > $archive/public-$today.all.masked.pb.gz
 ln -f `pwd`/public-$today.metadata.tsv.gz $archive/
 gzip -c public-$today.all.masked.nextclade.pangolin.pb \
     > $archive/public-$today.all.masked.nextclade.pangolin.pb.gz
 if [ -s lineageToPublicName ]; then
     gzip -c lineageToPublicName > $archive/lineageToPublicName.tsv.gz