014bee401f2fe14c79caf74d42b5f4f6cf797dd3 angie Thu Jun 18 15:23:11 2026 -0700 Loosen filters to accommodate increasingly divergent sequences. Tweak shell variables usherDir and friends to make it easier to try out different build directories. diff --git src/hg/utils/otto/sarscov2phylo/updateCombinedTree.sh src/hg/utils/otto/sarscov2phylo/updateCombinedTree.sh index 3cfa1cbde34..b38f5387d14 100755 --- src/hg/utils/otto/sarscov2phylo/updateCombinedTree.sh +++ src/hg/utils/otto/sarscov2phylo/updateCombinedTree.sh @@ -47,43 +47,44 @@ # Instead of the cluster, use Cheng's blazingly fast new usher-sampled: gunzip -f new.masked.mpl.gz time $usher \ -T 50 -A -e 5 \ -i prevRenamed.pb.gz \ --diff new.masked.mpl \ --ref $ottoDir/NC_045512.2.fa \ -o merged.pb.gz \ --optimization_radius 0 --batch_size_per_process 10 \ > usher.addNew.log 2>usher-sampled.stderr pigz -f -p 8 new.masked.mpl # Branch-specific masking time $scriptDir/maskDelta.sh merged.pb.gz merged.deltaMasked.pb.gz # Prune samples with too many private mutations and internal branches that are too long. $matUtils extract -i merged.deltaMasked.pb.gz \ - --max-parsimony 20 \ - --max-branch-length 70 \ - --max-path-length 225 \ + --max-parsimony 25 \ + --max-branch-length 75 \ + --max-path-length 250 \ -O -o merged.deltaMasked.filtered.pb.gz # matOptimize: used -r 8 -M2 until 2023-05-12, then switched to Cheng's recommended # -m 0.00000001 -M 4 (avoid identical-child-node problem in # https://github.com/sars-cov-2-variants/lineage-proposals/issues/40) # The -M 4 allowed up to radius 32, and crazy things started happening all while I was # trying to get the tree cleaned up for pango-designation release 1.20 --> lineageTree. # After 2023-05-20, when I found that matOptimize had moved a big chunk of B.1 onto a # B.1.1.7 garbage branch, causing big trouble for lineageTree (23_05_18_updateLineageTreePb.txt). # After that I changed it back to -M 2 for my sanity. If the identical-child thing happens again, # then I'll probably just run matOptimize twice, with a small radius the second time. + cd $ottoDir/$today && $matOptimize \ -T 64 -m 0.00000001 -M 2 -S move_log.filtered \ -i merged.deltaMasked.filtered.pb.gz \ -o gisaidAndPublic.$today.masked.preTrim.pb.gz \ >& matOptimize.filtered.log # Fix grandparent-reversion nodes that cause some lineages to be incorrectly placed as # sublineages of siblings. $matUtils fix -i gisaidAndPublic.$today.masked.preTrim.pb.gz -c 10 \ -o gisaidAndPublic.$today.masked.preTrim.fix.pb.gz # Again prune samples with too many private mutations and internal branches that are too long. $matUtils extract -i gisaidAndPublic.$today.masked.preTrim.fix.pb.gz \ --max-parsimony 20 \ --max-branch-length 70 \ @@ -162,69 +163,69 @@ | sort > tmp2 paste <(zcat gisaidAndPublic.$today.metadata.tsv.gz | cut -f 1-9 | head -1) \ <(echo -e "Nextstrain_clade_usher\tpango_lineage_usher") \ > gisaidAndPublic.$today.metadata.tsv join -t$'\t' tmp1 tmp2 \ >> gisaidAndPublic.$today.metadata.tsv pigz -p 8 -f gisaidAndPublic.$today.metadata.tsv rm tmp1 tmp2 # EPI_ISL_ ID to public sequence name mapping, so if users upload EPI_ISL IDs for which we have # public names & IDs, we can match them. cut -f 1,3 $epiToPublic > epiToPublic.latest # Update links to latest public+GISAID protobuf and metadata in /gbdb/wuhCor1/hgPhyloPlaceData/ dir=/gbdb/wuhCor1/hgPhyloPlaceData -ln -sf `pwd`/gisaidAndPublic.$today.masked.pb.gz $dir/public.plusGisaid.latest.masked.pb.gz -ln -sf `pwd`/gisaidAndPublic.$today.metadata.tsv.gz $dir/public.plusGisaid.latest.metadata.tsv.gz -ln -sf `pwd`/hgPhyloPlace.plusGisaid.description.txt $dir/public.plusGisaid.latest.version.txt -ln -sf `pwd`/epiToPublic.latest $dir/ -ln -sf `pwd`/samples.$today.gz $dir/public.plusGisaid.names.gz +ssh hgwdev ln -sf `pwd`/gisaidAndPublic.$today.masked.pb.gz $dir/public.plusGisaid.latest.masked.pb.gz +ssh hgwdev ln -sf `pwd`/gisaidAndPublic.$today.metadata.tsv.gz $dir/public.plusGisaid.latest.metadata.tsv.gz +ssh hgwdev ln -sf `pwd`/hgPhyloPlace.plusGisaid.description.txt $dir/public.plusGisaid.latest.version.txt +ssh hgwdev ln -sf `pwd`/epiToPublic.latest $dir/ +ssh hgwdev ln -sf `pwd`/samples.$today.gz $dir/public.plusGisaid.names.gz # Extract Omicron-only tree for faster searches $matUtils extract -i gisaidAndPublic.$today.masked.pb.gz -c B.1.1.529 \ -o gisaidAndPublic.$today.masked.omicron.pb.gz \ -u samples.$today.omicron pigz -f -p 8 samples.$today.omicron -ln -sf $(pwd)/gisaidAndPublic.$today.masked.omicron.pb.gz \ +ssh hgwdev ln -sf $(pwd)/gisaidAndPublic.$today.masked.omicron.pb.gz \ /gbdb/wuhCor1/hgPhyloPlaceData/public.plusGisaid.latest.masked.omicron.pb.gz sampleCountComma=$(echo $(zcat samples.$today.omicron.gz | wc -l) \ | sed -re 's/([0-9]+)([0-9]{3})$/\1,\2/; s/([0-9]+)([0-9]{3},[0-9]{3})$/\1,\2/;') echo "$sampleCountComma Omicron genomes from GISAID, GenBank, COG-UK and CNCB ($today); sarscov2phylo 13-11-20 tree with newer sequences added by UShER" \ > hgPhyloPlace.plusGisaid.omicron.description.txt -ln -sf $(pwd)/hgPhyloPlace.plusGisaid.omicron.description.txt \ +ssh hgwdev ln -sf $(pwd)/hgPhyloPlace.plusGisaid.omicron.description.txt \ /gbdb/wuhCor1/hgPhyloPlaceData/public.plusGisaid.latest.omicron.version.txt # Memory-mapped hash tables for metadata and name lookup tabToMmHash gisaidAndPublic.$today.metadata.tsv.gz gisaidAndPublic.$today.metadata.mmh -ln -sf $(pwd)/gisaidAndPublic.$today.metadata.mmh \ +ssh hgwdev ln -sf $(pwd)/gisaidAndPublic.$today.metadata.mmh \ /gbdb/wuhCor1/hgPhyloPlaceData/public.plusGisaid.latest.metadata.mmh for s in samples.$today.gz samples.$today.omicron.gz; do zcat $s \ | awk -F\| '{ print $0 "\t" $0; print $1 "\t" $0; if ($3 != "") { print $2 "\t" $0; } }' \ | tawk '$1 != "RNA" && $1 !~ /\/RNA\// && $1 !~/^Germany\/Molecular_surveillance_of_SARS/ && \ $1 !~ /^Iceland\/SARS-CoV-2_Iceland/' \ > nameLookup.tab cut -f 1,3 $epiToPublic \ | subColumn -skipMiss 2 stdin nameLookup.tab tmp.tab cat tmp.tab >> nameLookup.tab rm tmp.tab mmh=$(basename $s .gz).mmh tabToMmHash nameLookup.tab $mmh rm nameLookup.tab pub=$(echo $mmh | sed -re 's/'samples.$today'/public.plusGisaid.names/') - ln -sf $(pwd)/$mmh /gbdb/wuhCor1/hgPhyloPlaceData/$pub + ssh hgwdev ln -sf $(pwd)/$mmh /gbdb/wuhCor1/hgPhyloPlaceData/$pub done # Make Taxonium v2 protobuf for display usher_to_taxonium --input gisaidAndPublic.$today.masked.pb.gz \ --metadata gisaidAndPublic.$today.metadata.tsv.gz \ --genbank ~angie/github/taxonium/taxoniumtools/test_data/hu1.gb \ --columns genbank_accession,country,date,pangolin_lineage,pango_lineage_usher \ --clade_types=nextstrain,pango \ --name_internal_nodes \ --title "$today tree with sequences from GISAID, INSDC, COG-UK and CNCB" \ --output gisaidAndPublic.$today.masked.taxonium.jsonl.gz \ >& utt.log & $scriptDir/extractPublicTree.sh $today $prevDate