8d5f0115cf9ab2b99f18a1467c0be02a6ae55de4 angie Thu Jun 18 15:18:03 2026 -0700 Use taxonium overlay html to describe the trees. Adjust filters. 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/fluA/buildConcatTreeD1.1.sh src/hg/utils/otto/fluA/buildConcatTreeD1.1.sh index 184fe05f541..af278e77860 100755 --- src/hg/utils/otto/fluA/buildConcatTreeD1.1.sh +++ src/hg/utils/otto/fluA/buildConcatTreeD1.1.sh @@ -141,31 +141,33 @@ cladeNode=$(grep $(segRefAcc $seg) sample-paths.$asmAcc.$segRef \ | awk '{print $NF;}' \ | sed -re 's/:.*//;') grep -w $cladeNode sample-paths.$asmAcc.$segRef \ | cut -f 1 > samples.h5n1_D1.1_2024.$seg # Use GenBank sequences found in the trees plus Andersen Lab assembled sequences. cat samples.h5n1_D1.1_2024.* | grep -v \|SRR | cut -d\| -f 2 \ | grep -Fwf - $fluANcbiDir/metadata.tsv \ | grep -v 1969-12-31 \ | grep -v 1970-01-01 \ | cut -f 2,15 \ | tawk '{ if ($2 != "") { print $2; } else { print $1; } }' \ | sort -u \ | grep -Ff - $fluANcbiDir/metadata.tsv \ +| grep -vE ' clone ?[0-9]+' \ | cut -f 1,17 \ +| grep -Fwf <(zcat renaming.tsv.gz | cut -f 1) \ > cladeAccToSeg # Extract the sequences into per-segment fasta files... renamed from accession to tree name. # joinSegments.py below will ignore the uniquifying INSDC accession part of names. Remove the # uniquifying segment name from Andersen Lab sequences. for seg in 1 2 3 4 5 6 7 8; do tawk '$2 == '$seg' {print $1;}' cladeAccToSeg \ | faSomeRecords <(xzcat $fluANcbiDir/genbank.fa.xz) stdin stdout \ | faRenameRecords stdin renaming.tsv.gz h5n1_D1.1_2024.$seg.fa segName=$(segName $seg) set +o pipefail fastaNames $fluADir/andersen_lab.srrNotGb.renamed.fa \ | grep _$segName \ | grep -Fwf <(grep \|SRR samples.h5n1_D1.1_2024.* | cut -d\| -f 2) \ | faSomeRecords $fluADir/andersen_lab.srrNotGb.renamed.fa stdin stdout \ | sed -re '/^>/ s@_'$segName'/@/@;' \ @@ -260,61 +262,62 @@ haFields=$(zcat H5N1_HA_DMS_metadata.tsv.gz | cut -f 2- | head -1 | sed -re 's/\t/\\t/g') pb2Fields=$(zcat PB2_DMS_metadata.tsv.gz | cut -f 2-3 | head -1 | sed -re 's/\t/\\t/g') set -o pipefail echo -e "$oldFields\t$haFields\t$pb2Fields" > h5n1_D1.1_2024.metadata.tsv cat tmp >> h5n1_D1.1_2024.metadata.tsv rm tmp pigz -f -p 8 h5n1_D1.1_2024.metadata.tsv usher_to_taxonium --input h5n1_D1.1_2024.pb.gz \ --metadata h5n1_D1.1_2024.metadata.tsv.gz \ --columns host,country,location,date,authors,mouse_escape,ferret_escape,cell_entry,stability,sa26_increase,mouse_escape_mutations,ferret_escape_mutations,cell_entry_mutations,stability_mutations,sa26_increase_mutations,mutdiffsel,mutdiffsel_mutations \ --genbank $fluADir/h5n1_D1.1_2024/concat.gbff \ --name_internal_nodes \ --title "2024 H5N1 D1.1 outbreak in USA, concatenated segments from INSDC and SRA ($today)" \ --config_json $fluAScriptDir/concat.config.json \ + --overlay_html $fluAScriptDir/taxonium_overlay_d1_1.html \ --chronumental \ --chronumental_steps 500 \ --chronumental_add_inferred_date chronumental_date \ --output h5n1_D1.1_2024.jsonl.gz \ >& utt.log # Link to /gbdb/ location dir=/gbdb/wuhCor1/hgPhyloPlaceData/influenzaA/h5n1_D1.1_2024 -mkdir -p $dir -ln -sf $(pwd)/h5n1_D1.1_2024.pb.gz $dir/h5n1_D1.1_2024.latest.pb.gz -ln -sf $(pwd)/h5n1_D1.1_2024.metadata.tsv.gz $dir/h5n1_D1.1_2024.latest.metadata.tsv.gz -ln -sf $(pwd)/hgPhyloPlace.description.h5n1_D1.1_2024.txt \ +ssh hgwdev mkdir -p $dir +ssh hgwdev ln -sf $(pwd)/h5n1_D1.1_2024.pb.gz $dir/h5n1_D1.1_2024.latest.pb.gz +ssh hgwdev ln -sf $(pwd)/h5n1_D1.1_2024.metadata.tsv.gz $dir/h5n1_D1.1_2024.latest.metadata.tsv.gz +ssh hgwdev ln -sf $(pwd)/hgPhyloPlace.description.h5n1_D1.1_2024.txt \ $dir/h5n1_D1.1_2024.latest.version.txt # Extract Newick and VCF for anyone who wants to download those instead of protobuf $matUtils extract -i h5n1_D1.1_2024.pb.gz \ -t h5n1_D1.1_2024.nwk \ -v h5n1_D1.1_2024.vcf >& tmp.log pigz -p 8 -f h5n1_D1.1_2024.nwk h5n1_D1.1_2024.vcf # Make a ref + all fasta download file for Delphy folks cat $fluADir/h5n1_D1.1_2024/concat.fa h5n1_D1.1_2024.aligned.fa \ | pigz -p 8 \ > h5n1_D1.1_2024.msa.fa.gz # Update 'latest' in $archiveRoot for e in jsonl.gz metadata.tsv.gz nwk.gz pb.gz vcf.gz msa.fa.gz ; do ln -sf $(pwd)/h5n1_D1.1_2024.$e $archiveRoot/h5n1_D1.1_2024.latest.$e done ln -sf $(pwd)/hgPhyloPlace.description.h5n1_D1.1_2024.txt \ $archiveRoot/h5n1_D1.1_2024.latest.version.txt # Update hgdownload-test link for archive (adding assembly/segRef hierarchy) -mkdir -p $downloadsRoot/$asmDir/UShER_h5n1_D1.1_2024 -ln -sf $archiveRoot/h5n1_D1.1_2024.latest.* $downloadsRoot/$asmDir/UShER_h5n1_D1.1_2024/ +ssh hgwdev mkdir -p $downloadsRoot/$asmDir/UShER_h5n1_D1.1_2024 +ssh hgwdev ln -sf $archiveRoot/h5n1_D1.1_2024.latest.* $downloadsRoot/$asmDir/UShER_h5n1_D1.1_2024/ # rsync to hgdownload hubs dir -for h in hgdownload1 hgdownload3; do - if rsync -a -L --delete $downloadsRoot/$asmDir/UShER_h5n1_D1.1_2024 \ +for h in hgdownload1 hgdownload2 hgdownload3; do + if ssh hgwdev rsync -a -L --delete $downloadsRoot/$asmDir/UShER_h5n1_D1.1_2024 \ qateam@$h:/mirrordata/hubs/$asmDir/ ; then true else echo "" echo "*** rsync to $h failed -- disk full ? ***" echo "" fi done