de23dd389f8e8e35c41880a85d74e9838a4df8e4 angie Thu Jun 18 15:11:49 2026 -0700 ssh to hgwdev for commands that need access to hgwdev-only directories, because these scripts run elsewhere now. Tweak shell variables usherDir and friends to make it easier to try out different build directories. diff --git src/hg/utils/otto/mpxv/buildTree.sh src/hg/utils/otto/mpxv/buildTree.sh index aafd231ec2b..fab2118ce73 100755 --- src/hg/utils/otto/mpxv/buildTree.sh +++ src/hg/utils/otto/mpxv/buildTree.sh @@ -1,255 +1,254 @@ #!/bin/bash source ~/.bashrc set -beEu -o pipefail mpxvScriptDir=$(dirname "${BASH_SOURCE[0]}") today=$(date +%F) mpxvDir=/hive/data/outside/otto/mpxv # RefSeq assembly for our track hubs for clade I and clade II asmAcc_I=GCF_000857045.1 asmAcc_II=GCF_014621545.1 # GenBank flat file for taxonium gbff_I=$mpxvDir/NC_003310.1.gbff gbff_II=$mpxvDir/NC_063383.1.gbff # Metadata file metadata_I=$mpxvNcbiDir/metadata.cladeI.tsv metadata_II=$mpxvNcbiDir/metadata.2017outbreak.tsv # Mask file #*** TODO get mask regions for clade I! mask_I=/dev/null mask_II=$mpxvScriptDir/mask.vcf.gz # faToVcf -maxDiff maxDiff_I=10000 maxDiff_II=100 # Minimum number of samples in tree to be sure the NCBI download wasn't incomplete minSamples_I=600 minSamples_II=7000 # Download archive storage directory archiveRoot_I=/hive/users/angie/publicTreesMpxvCladeI archiveRoot_II=/hive/users/angie/publicTreesHMPXV # Download archive URL directory downloadDir_I=UShER_MPXV_cladeI downloadDir_II=UShER_hMPXV mpxvNcbiDir=$mpxvDir/ncbi/ncbi.latest -usherDir=~angie/github/usher -usherSampled=$usherDir/build/usher-sampled -usher=$usherDir/build/usher -matUtils=$usherDir/build/matUtils +usherDir=~angie/github/usher/build +usherSampled=$usherDir/usher-sampled +matUtils=$usherDir/matUtils +matOptimize=$usherDir/matOptimize if [[ ! -d $mpxvDir/ncbi/ncbi.$today || ! -s $mpxvDir/ncbi/ncbi.$today/genbank.fa.xz ]]; then mkdir -p $mpxvDir/ncbi/ncbi.$today $mpxvScriptDir/getNcbiMpxv.sh >& $mpxvDir/ncbi/ncbi.$today/getNcbiMpxv.log fi buildDir=$mpxvDir/build/$today mkdir -p $buildDir cd $buildDir # This builds the whole tree from scratch! Eventually we'll want to add only the new sequences # to yesterday's tree. echo '()' > emptyTree.nwk for clade in I II; do if [[ $clade == "I" ]]; then asmAcc=$asmAcc_I gbff=$gbff_I metadata=$metadata_I mask=$mask_I maxDiff=$maxDiff_I minSamples=$minSamples_I archiveRoot=$archiveRoot_I downloadDir=$downloadDir_I else asmAcc=$asmAcc_II gbff=$gbff_II metadata=$metadata_II mask=$mask_II maxDiff=$maxDiff_II minSamples=$minSamples_II archiveRoot=$archiveRoot_II downloadDir=$downloadDir_II fi # Use metadata to make a renaming file perl -wne 'chomp; @w=split(/\t/); my ($acc, $iso, $loc, $date, $name) = ($w[0], $w[1], $w[3], $w[4], $w[11]); if ($iso eq "") { $iso = $name; } my $country = $loc; $country =~ s/:.*//; my $COU = $country; $COU =~ s/^(\w{3}).*/$1/; $COU = uc($COU); if ($country eq "United Kingdom") { $COU = "UK"; } if ($iso !~ /$country/ && $iso !~ /\b$COU\b/) { $iso = "$country/$iso"; } my $year = $date; $year =~ s/-.*//; if ($iso !~ /$year/) { $iso = "$iso/$year"; } my $fullName = $name ? "$name|$acc|$date" : "$acc|$date"; $fullName =~ s/[ ,:()]/_/g; print "$acc\t$fullName\n";' \ $metadata > renaming.$clade.tsv # Use Nextstrain's exclusion list too. awk '{print $1;}' ~angie/github/monkeypox/phylogenetic/defaults/exclude_accessions.txt \ | grep -Fwf - <(cut -f 1 $metadata) \ > exclude.ids asmDir=$(echo $asmAcc \ | sed -re 's@^(GC[AF])_([0-9]{3})([0-9]{3})([0-9]{3})\.([0-9]+)@\1/\2/\3/\4/\1_\2\3\4.\5@') ref2bit=/hive/data/genomes/asmHubs/$asmDir/$asmAcc.2bit time cat <(twoBitToFa $ref2bit stdout) <(xzcat $mpxvNcbiDir/nextalign.$clade.fa.xz) \ | faToVcf -maxDiff=$maxDiff -verbose=2 -excludeFile=exclude.ids -includeNoAltN stdin stdout \ | vcfRenameAndPrune stdin renaming.$clade.tsv stdout \ | vcfFilter -excludeVcf=$mask stdin \ | pigz -p 8 \ > all.masked.vcf.gz time $usherSampled -T 16 -A -e 10 \ -t emptyTree.nwk \ -v all.masked.vcf.gz \ -o mpxv.clade$clade.$today.masked.preOpt.pb.gz \ --optimization_radius 0 --batch_size_per_process 10 \ > usher.addNew.log 2>usher-sampled.stderr # Optimize: - time ~angie/github/usher_branch/build/matOptimize \ - -T 16 -r 20 -M 2 -S move_log.usher_branch \ + time $matOptimize \ + -T 16 -r 20 -M 2 -S move_log \ -i mpxv.clade$clade.$today.masked.preOpt.pb.gz \ -o mpxv.clade$clade.$today.masked.opt.pb.gz \ - >& matOptimize.usher_branch.log + >& matOptimize.log # It crashes when I add # -v all.masked.vcf.gz \ # -- bug Cheng later. if [[ $clade == "II" ]]; then # Annotate root nodes for Nextstrain lineages. join -t$'\t' <(sort renaming.$clade.tsv) <(zcat $mpxvNcbiDir/nextclade.$clade.tsv.gz | cut -f 2,5 | sort) \ | tawk '{print $3, $2;}' | sort > lineageToName $matUtils annotate -T 16 -i mpxv.clade$clade.$today.masked.opt.pb.gz -c lineageToName \ -o mpxv.clade$clade.$today.masked.pb.gz \ >& annotate.$clade.log # Make metadata that uses same names as tree and includes nextclade lineage assignments. echo -e "strain\tgenbank_accession\tdate\tcountry\tlocation\tlength\thost\tbioproject_accession\tbiosample_accession\tsra_accession\tauthors\tpublications\tNextstrain_lineage" \ > mpxv.clade$clade.$today.metadata.tsv join -t$'\t' <(sort renaming.$clade.tsv) <(zcat $mpxvNcbiDir/nextclade.$clade.tsv.gz | cut -f 2,5 | sort) \ | join -t$'\t' -o 1.2,2.1,2.6,2.4,2.5,2.8,2.9,2.10,2.11,2.13,2.14,2.15,1.3 \ - <(sort $metadata \ | perl -F'/\t/' -walne '$F[3] =~ s/(: ?|$)/\t/; print join("\t", @F);') \ >> mpxv.clade$clade.$today.metadata.tsv else # No lineages to annotate; just clean up .opt with -O. $matUtils extract -i mpxv.clade$clade.$today.masked.opt.pb.gz \ -O -o mpxv.clade$clade.$today.masked.pb.gz # Make metadata that uses same names as tree and includes nextclade clade Ia/Ib assignments. echo -e "strain\tgenbank_accession\tdate\tcountry\tlocation\tlength\thost\tbioproject_accession\tbiosample_accession\tsra_accession\tauthors\tpublications\tNextstrain_clade" \ > mpxv.clade$clade.$today.metadata.tsv join -t$'\t' <(sort renaming.$clade.tsv) <(zcat $mpxvNcbiDir/nextclade.$clade.tsv.gz | cut -f 2,3 | sort) \ | join -t$'\t' -o 1.2,2.1,2.6,2.4,2.5,2.8,2.9,2.10,2.11,2.13,2.14,2.15,1.3 \ - <(sort $metadata \ | perl -F'/\t/' -walne '$F[3] =~ s/(: ?|$)/\t/; print join("\t", @F);') \ >> mpxv.clade$clade.$today.metadata.tsv fi pigz -f -p 8 mpxv.clade$clade.$today.metadata.tsv # Make a tree version description for hgPhyloPlace $matUtils extract -i mpxv.clade$clade.$today.masked.pb.gz -u samples.$clade.$today sampleCount=$(wc -l < samples.$clade.$today) # Sometimes NCBI download is incomplete; don't replace yesterday's tree in that case. if (( $sampleCount < $minSamples )); then echo "*** Too few samples ($sampleCount) for clade $clade! Expected at least $minSamples. Halting. ***" exit 1 fi sampleCountComma=$(echo $sampleCount \ | 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 INSDC (GenBank/ENA/DDBJ) ($today)" \ > hgPhyloPlace.description.$clade.txt # Make a taxonium view if [[ $clade == "I" ]]; then columns=genbank_accession,location,date,authors,Nextstrain_clade else columns=genbank_accession,location,date,authors,Nextstrain_lineage fi usher_to_taxonium --input mpxv.clade$clade.$today.masked.pb.gz \ --metadata mpxv.clade$clade.$today.metadata.tsv.gz \ --genbank $gbff \ --columns $columns \ --clade_types=pango \ --output mpxv.clade$clade.$today.masked.taxonium.jsonl.gz \ >& usher_to_taxonium.log # Update links to latest protobuf and metadata in /gbdb directories nc=$(basename $gbff .gbff) dir=/gbdb/wuhCor1/hgPhyloPlaceData/mpxv/$nc - mkdir -p $dir - ln -sf $(pwd)/mpxv.clade$clade.$today.masked.pb.gz $dir/mpxv.clade$clade.latest.pb.gz - ln -sf $(pwd)/mpxv.clade$clade.$today.metadata.tsv.gz $dir/mpxv.clade$clade.latest.metadata.tsv.gz - ln -sf $(pwd)/hgPhyloPlace.description.$clade.txt $dir/mpxv.clade$clade.latest.version.txt + ssh hgwdev mkdir -p $dir + ssh hgwdev ln -sf $(pwd)/mpxv.clade$clade.$today.masked.pb.gz $dir/mpxv.clade$clade.latest.pb.gz + ssh hgwdev ln -sf $(pwd)/mpxv.clade$clade.$today.metadata.tsv.gz $dir/mpxv.clade$clade.latest.metadata.tsv.gz + ssh hgwdev ln -sf $(pwd)/hgPhyloPlace.description.$clade.txt $dir/mpxv.clade$clade.latest.version.txt # Extract Newick and VCF for anyone who wants to download those instead of protobuf $matUtils extract -i mpxv.clade$clade.$today.masked.pb.gz \ -t mpxv.clade$clade.$today.nwk \ -v mpxv.clade$clade.$today.masked.vcf pigz -p 8 -f mpxv.clade$clade.$today.nwk mpxv.clade$clade.$today.masked.vcf # Link to public trees download directory hierarchy read y m d < <(echo $today | sed -re 's/-/ /g') archive=$archiveRoot/$y/$m/$d mkdir -p $archive if [[ $clade == "I" ]]; then ln -f $(pwd)/mpxv.clade$clade.$today.{nwk,masked.vcf,metadata.tsv,masked.taxonium.jsonl,masked.pb}.gz $archive/ ln -f $(pwd)/hgPhyloPlace.description.$clade.txt $archive/mpxv.clade$clade.$today.version.txt else # At first, clade II hMPXV was the only mpox, so cladeII was not part of the download file names. # Also, for clade II we're not making a tree with all clade II, only the 2017 outbreak, so cladeII # in the name wouldn't be entirely accurate either. for f in mpxv.clade$clade.$today.{nwk,masked.vcf,metadata.tsv,masked.taxonium.jsonl,masked.pb}.gz; do downloadF=$(echo $f | sed -re 's/.cladeII//;') ln -f $(pwd)/$f $archive/$downloadF done ln -f $(pwd)/hgPhyloPlace.description.$clade.txt $archive/mpxv.$today.version.txt fi # Update 'latest' in $archiveRoot for f in $archive/mpxv*.$today.*; do latestF=$(echo $(basename $f) | sed -re 's/'$today'/latest/') ln -f $f $archiveRoot/$latestF done # Update hgdownload-test link for archive - mkdir -p /usr/local/apache/htdocs-hgdownload/hubs/$asmDir/$downloadDir/$y/$m - ln -sf $archive /usr/local/apache/htdocs-hgdownload/hubs/$asmDir/$downloadDir/$y/$m + ssh hgwdev ln -sf $archiveRoot/*.latest.* /data/apache/htdocs-hgdownload/hubs/$asmDir/$downloadDir/ # rsync to hgdownload hubs dir - for h in hgdownload1 hgdownload3; do - if rsync -a -L --delete /usr/local/apache/htdocs-hgdownload/hubs/$asmDir/$downloadDir/* \ + for h in hgdownload1 hgdownload2 hgdownload3; do + if ssh hgwdev rsync -a -L --delete /data/apache/htdocs-hgdownload/hubs/$asmDir/$downloadDir/* \ qateam@$h:/mirrordata/hubs/$asmDir/$downloadDir/; then true else echo "" echo "*** rsync to $h failed; disk full? ***" echo "" fi done if [[ $clade == "II" ]]; then set +o pipefail grep 'Could not' annotate.$clade.log | cat grep skipping annotate.$clade.log | cat set -o pipefail fi cat hgPhyloPlace.description.$clade.txt zcat mpxv.clade$clade.$today.metadata.tsv.gz | tail -n+2 | cut -f 13 | sort | uniq -c done