6091af3e8fff91595ab6361974787afef032a342 angie Thu Jun 25 11:03:38 2026 -0700 Fix second filtering to match first. Extract tree of sequences from the past year plus sequences within a radius of 5 edges of those. diff --git src/hg/utils/otto/sarscov2phylo/updateCombinedTree.sh src/hg/utils/otto/sarscov2phylo/updateCombinedTree.sh index b38f5387d14..61bb0df9542 100755 --- src/hg/utils/otto/sarscov2phylo/updateCombinedTree.sh +++ src/hg/utils/otto/sarscov2phylo/updateCombinedTree.sh @@ -74,33 +74,33 @@ # 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 \ - --max-path-length 225 \ + --max-parsimony 25 \ + --max-branch-length 75 \ + --max-path-length 250 \ -O -o gisaidAndPublic.$today.masked.pb.gz fi # Exclude sequences with a very high number of EPPs from future runs grep ^Current usher.addNew.log \ | awk '$16 >= 10 {print $8;}' \ | awk -F\| '{ if ($3 == "") { print $1; } else { print $2; } }' \ > tooManyEpps.ids cat tooManyEpps.ids >> ../tooManyEpps.ids $matUtils extract -i gisaidAndPublic.$today.masked.pb.gz -u samples.$today $scriptDir/combineMetadata.sh $prevDate $today # version/description files @@ -183,35 +183,55 @@ # 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 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 ssh hgwdev ln -sf $(pwd)/hgPhyloPlace.plusGisaid.omicron.description.txt \ /gbdb/wuhCor1/hgPhyloPlaceData/public.plusGisaid.latest.omicron.version.txt +# Extract recent-only tree (samples from the past year plus samples within radius 5 of those) for faster searches +set +x +conda activate bte +set -x +python $scriptDir/find_recent_samples.py -i gisaidAndPublic.$today.metadata.tsv.gz -n 365 -o recent_samples.txt +python $scriptDir/prune_to_radius.py -i gisaidAndPublic.$today.masked.pb.gz -l recent_samples.txt -r 5 \ + -o samples.$today.recent +$matUtils extract -i gisaidAndPublic.$today.masked.pb.gz \ + -s samples.$today.recent \ + -o gisaidAndPublic.$today.masked.recent.pb.gz +pigz -f -p 8 samples.$today.recent +ssh hgwdev ln -sf $(pwd)/gisaidAndPublic.$today.masked.recent.pb.gz \ + /gbdb/wuhCor1/hgPhyloPlaceData/public.plusGisaid.latest.masked.recent.pb.gz +sampleCountComma=$(echo $(zcat samples.$today.recent.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 recent genomes from GISAID, GenBank, COG-UK and CNCB ($today); sarscov2phylo 13-11-20 tree with newer sequences added by UShER" \ + > hgPhyloPlace.plusGisaid.recent.description.txt +ssh hgwdev ln -sf $(pwd)/hgPhyloPlace.plusGisaid.recent.description.txt \ + /gbdb/wuhCor1/hgPhyloPlaceData/public.plusGisaid.latest.recent.version.txt + # Memory-mapped hash tables for metadata and name lookup tabToMmHash gisaidAndPublic.$today.metadata.tsv.gz 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 +for s in samples.$today.gz samples.$today.omicron.gz samples.$today.recent.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/') ssh hgwdev ln -sf $(pwd)/$mmh /gbdb/wuhCor1/hgPhyloPlaceData/$pub done