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/buildTree.sh src/hg/utils/otto/fluA/buildTree.sh index cf5130f2d8e..61420ef5aa0 100755 --- src/hg/utils/otto/fluA/buildTree.sh +++ src/hg/utils/otto/fluA/buildTree.sh @@ -4,35 +4,34 @@ # Align INSDC sequences to reference and build N*M trees where N = 8 (number of segments in the # influenza genome) and M = 7 (number of RefSeq genome assemblies) fluAScriptDir=$(dirname "${BASH_SOURCE[0]}") if [[ $# > 0 ]]; then today=$1 else today=$(date +%F) fi fluADir=/hive/data/outside/otto/fluA fluANcbiDir=$fluADir/ncbi/ncbi.latest -usherDir=~angie/github/usher -usherSampled=$usherDir/build/usher-sampled -usher=$usherDir/build/usher -matUtils=$usherDir/build/matUtils -matOptimize=$usherDir/build/matOptimize +usherDir=~angie/github/usher/build +usherSampled=$usherDir/usher-sampled +matUtils=$usherDir/matUtils +matOptimize=$usherDir/matOptimize minSize=800 threads=16 assemblyDir=/hive/data/outside/ncbi/genomes asmHubDir=/hive/data/genomes/asmHubs/refseqBuild archiveRoot=/hive/users/angie/publicTreesFluA downloadsRoot=/data/apache/htdocs-hgdownload/hubs # assembly serotype taxid isolate # GCF_000865085.1 H3N2 335341 A/New York/392/2004(H3N2) # GCF_001343785.1 H1N1 641809 A/California/07/2009(H1N1) # GCF_000865725.1 H1N1 211044 A/Puerto Rico/8/1934(H1N1) # GCF_000928555.1 H7N9 1332244 A/Shanghai/02/2013(H7N9) @@ -239,56 +238,59 @@ >& tmp.log awk -F\| '{if ($3 == "") { print $1; } else { print $2; }}' samples.$asmAcc.$segRef.$today \ > accs.$asmAcc.$segRef.tsv sampleCountComma=$(wc -l < samples.$asmAcc.$segRef.$today \ | sed -re 's/([0-9]+)([0-9]{3})$/\1,\2/; s/([0-9]+)([0-9]{3},[0-9]{3})$/\1,\2/;') if [[ $asmAcc == GCF_000864105.1 ]]; then echo "$sampleCountComma genomes from INSDC (GenBank/ENA/DDBJ) and/or GISAID ($today)" \ > hgPhyloPlace.description.$asmAcc.$segRef.txt else echo "$sampleCountComma genomes from INSDC (GenBank/ENA/DDBJ) ($today)" \ > hgPhyloPlace.description.$asmAcc.$segRef.txt fi # Depending on the segment RefSeq, maybe run nextclade - # Note: nextclade has a dataset flu_h3n2_na but it does not assign clades. case $segRef in "NC_007366.1") nextcladeName=flu_h3n2_ha ;; + "NC_007368.1") + nextcladeName=flu_h3n2_na + ;; "NC_026433.1") nextcladeName=flu_h1n1pdm_ha ;; "NC_026434.1") nextcladeName=flu_h1n1pdm_na ;; "NC_007362.1") nextcladeName=community/moncla-lab/iav-h5/ha/all-clades ;; *) nextcladeName="" nextcladeTaxCo="" ;; esac if [[ x$nextcladeName != x ]]; then nextclade dataset get --name $nextcladeName --output-zip $nextcladeName.zip (if [[ $segRef == "NC_007362.1" ]]; then # Also run on collab's sequences and SRA assemblies time cat <(faSomeRecords <(xzcat $fluADir/ncbi/ncbi.$today/genbank.fa.xz) \ accs.$asmAcc.$segRef.tsv stdout) \ - $fluADir/h5nx.epiNoMatchRenamed.fa \ $fluADir/andersen_lab.srrNotGb.renamed.fa + $fluADir/h5nx.epiNoMatchRenamed.fa \ + $fluADir/andersen_lab.srrNotGb.renamed.fa else time faSomeRecords <(xzcat $fluADir/ncbi/ncbi.$today/genbank.fa.xz) \ accs.$asmAcc.$segRef.tsv stdout fi) \ | nextclade run \ -D $nextcladeName.zip \ -j $threads \ --retry-reverse-complement true \ --output-tsv nextclade.$asmAcc.$segRef.tsv \ --output-columns-selection seqName,clade,totalSubstitutions,totalDeletions,totalInsertions,totalMissing,totalNonACGTNs,alignmentStart,alignmentEnd,substitutions,deletions,insertions,aaSubstitutions,aaDeletions,aaInsertions,missing,unknownAaRanges,nonACGTNs \ >& nextclade.$asmAcc.$segRef.log nextcladeTaxCo=",Nextstrain_clade" fi # Make metadata that uses same names as tree @@ -398,44 +400,54 @@ - <(tail -n+2 PB2_DMS_metadata.tsv | sort) \ >> fluA.$asmAcc.$segRef.$today.metadata.tsv fi wc -l fluA.$asmAcc.$segRef.$today.metadata.tsv pigz -f -p 8 fluA.$asmAcc.$segRef.$today.metadata.tsv # Make a taxonium view if [[ $asmAcc == GCF_000864105.1 ]]; then gisaid=" and/or GISAID" else gisaid="" fi title="Influenza A $strain segment $segment ($segName) $today tree with $sampleCountComma genomes from INSDC$gisaid" columns="genbank_accession,country,location,date,host,serotype,segment,genoflu_genotype,genoflu_segtype,authors,bioproject_accession$nextcladeTaxCo" config_json="" + overlay_template=$fluAScriptDir/taxonium_overlay.html if [[ $segRef == "NC_007362.1" ]]; then columns="$columns,mouse_escape,ferret_escape,cell_entry,stability,sa26_increase" config_json="--config_json $fluAScriptDir/h5n1_ha.config.json" + overlay_template=$fluAScriptDir/taxonium_overlay_h5n1_ha.html elif [[ $segment == 1 ]]; then columns="$columns,mutdiffsel,mutdiffsel_mutations,aa_substitution_count" config_json="--config_json $fluAScriptDir/pb2.config.json" + overlay_template=$fluAScriptDir/taxonium_overlay_pb2.html fi + sed -re 's@__ASMDIR__@'$asmDir'@g' $overlay_template \ + | sed -re 's@__STRAIN__@'"$strain"'@g' \ + | sed -re 's/__SEGNUM__/'$segment'/g' \ + | sed -re 's/__SEGNAME__/'$segName'/g' \ + | sed -re 's/__SEGREF__/'$segRef'/g' \ + > overlay.html if \ usher_to_taxonium --input fluA.$asmAcc.$segRef.$today.pb \ --metadata fluA.$asmAcc.$segRef.$today.metadata.tsv.gz \ --columns $columns \ --genbank $fluADir/$asmAcc/$segRef.gbff \ --name_internal_nodes \ $config_json \ + --overlay_html overlay.html \ --title "$title" \ --output fluA.$asmAcc.$segRef.$today.taxonium.jsonl.gz \ >& utt.log; then true; else mv utt.log utt.$asmAcc.$segRef.fail.log echo "*** usher_to_taxonium failed, see utt.$asmAcc.$segRef.fail.log ***" fi pigz -f -p 8 samples.$asmAcc.$segRef.$today if [[ $asmAcc == GCF_000864105.1 ]]; then # Rename so that we omit collab-shared sequences from files named public, # and make truly public versions by pruning down the tree mv accs.$asmAcc.$segRef{,.plusGisaid}.tsv @@ -456,118 +468,126 @@ set +o pipefail zcat fluA.$asmAcc.$segRef.plusGisaid.$today.metadata.tsv.gz \ | head -1 \ > tmp set -o pipefail zcat fluA.$asmAcc.$segRef.plusGisaid.$today.metadata.tsv.gz \ | grep -Fwf accs.$asmAcc.$segRef.tsv \ >> tmp mv tmp fluA.$asmAcc.$segRef.$today.metadata.tsv pigz -f -p 8 fluA.$asmAcc.$segRef.$today.metadata.tsv sampleCountComma=$(zcat samples.$asmAcc.$segRef.$today.gz | wc -l \ | sed -re 's/([0-9]+)([0-9]{3})$/\1,\2/; s/([0-9]+)([0-9]{3},[0-9]{3})$/\1,\2/;') title="Influenza A $strain segment $segment ($segName) $today tree with $sampleCountComma genomes from INSDC" columns="genbank_accession,country,location,date,host,serotype,segment,genoflu_genotype,genoflu_segtype,authors,bioproject_accession$nextcladeTaxCo" config_json="" + overlay_template=$fluAScriptDir/taxonium_overlay.html if [[ $segRef == "NC_007362.1" ]]; then columns="$columns,mouse_escape,ferret_escape,cell_entry,stability,sa26_increase" config_json="--config_json $fluAScriptDir/h5n1_ha.config.json" + overlay_template=$fluAScriptDir/taxonium_overlay_h5n1_ha.html elif [[ $segRef == "NC_007357.1" ]]; then columns="$columns,mutdiffsel,mutdiffsel_mutations,aa_substitution_count" config_json="--config_json $fluAScriptDir/pb2.config.json" + overlay_template=$fluAScriptDir/taxonium_overlay_pb2.html fi + sed -re 's@__ASMDIR__@'$asmDir'@g' $overlay_template \ + | sed -re 's@__STRAIN__@'"$strain"'@g' \ + | sed -re 's/__SEGNUM__/'$segment'/g' \ + | sed -re 's/__SEGNAME__/'$segName'/g' \ + | sed -re 's/__SEGREF__/'$segRef'/g' \ + > overlay.html if \ usher_to_taxonium --input fluA.$asmAcc.$segRef.$today.pb \ --metadata fluA.$asmAcc.$segRef.$today.metadata.tsv.gz \ --columns $columns \ --genbank $fluADir/$asmAcc/$segRef.gbff \ --name_internal_nodes \ $config_json \ + --overlay_html overlay.html \ --title "$title" \ --output fluA.$asmAcc.$segRef.$today.taxonium.jsonl.gz \ >& utt.log; then true; else mv utt.log utt.$asmAcc.$segRef.fail.log echo "*** usher_to_taxonium failed, see utt.$asmAcc.$segRef.fail.log ***" fi echo "$sampleCountComma genomes from INSDC (GenBank/ENA/DDBJ) ($today)" \ > hgPhyloPlace.description.$asmAcc.$segRef.txt # Put .plusGisaid versions in non-hgdownload location dir=/gbdb/wuhCor1/hgPhyloPlaceData/influenzaA/$segRef - mkdir -p $dir - ln -sf $(pwd)/fluA.$asmAcc.$segRef.plusGisaid.$today.pb \ + ssh hgwdev mkdir -p $dir + ssh hgwdev ln -sf $(pwd)/fluA.$asmAcc.$segRef.plusGisaid.$today.pb \ $dir/fluA.$asmAcc.$segRef.plusGisaid.latest.pb - ln -sf $(pwd)/fluA.$asmAcc.$segRef.plusGisaid.$today.metadata.tsv.gz \ + ssh hgwdev ln -sf $(pwd)/fluA.$asmAcc.$segRef.plusGisaid.$today.metadata.tsv.gz \ $dir/fluA.$asmAcc.$segRef.plusGisaid.latest.metadata.tsv.gz - ln -sf $(pwd)/hgPhyloPlace.description.$asmAcc.$segRef.plusGisaid.txt \ + ssh hgwdev ln -sf $(pwd)/hgPhyloPlace.description.$asmAcc.$segRef.plusGisaid.txt \ $dir/fluA.$asmAcc.$segRef.plusGisaid.latest.version.txt - ln -sf $(pwd)/samples.$asmAcc.$segRef.plusGisaid.$today.gz \ + ssh hgwdev ln -sf $(pwd)/samples.$asmAcc.$segRef.plusGisaid.$today.gz \ $dir/fluA.$asmAcc.$segRef.plusGisaid.latest.samples.gz fi # Link regular versions to non-hgdownload location too dir=/gbdb/wuhCor1/hgPhyloPlaceData/influenzaA/$segRef - mkdir -p $dir - ln -sf $(pwd)/fluA.$asmAcc.$segRef.$today.pb \ + ssh hgwdev mkdir -p $dir + ssh hgwdev ln -sf $(pwd)/fluA.$asmAcc.$segRef.$today.pb \ $dir/fluA.$asmAcc.$segRef.latest.pb - ln -sf $(pwd)/fluA.$asmAcc.$segRef.$today.metadata.tsv.gz \ + ssh hgwdev ln -sf $(pwd)/fluA.$asmAcc.$segRef.$today.metadata.tsv.gz \ $dir/fluA.$asmAcc.$segRef.latest.metadata.tsv.gz - ln -sf $(pwd)/hgPhyloPlace.description.$asmAcc.$segRef.txt \ + ssh hgwdev ln -sf $(pwd)/hgPhyloPlace.description.$asmAcc.$segRef.txt \ $dir/fluA.$asmAcc.$segRef.latest.version.txt - ln -sf $(pwd)/samples.$asmAcc.$segRef.$today.gz \ + ssh hgwdev ln -sf $(pwd)/samples.$asmAcc.$segRef.$today.gz \ $dir/fluA.$asmAcc.$segRef.latest.samples.gz # Extract Newick and VCF for anyone who wants to download those instead of protobuf $matUtils extract -i fluA.$asmAcc.$segRef.$today.pb \ -t fluA.$asmAcc.$segRef.$today.nwk \ -v fluA.$asmAcc.$segRef.$today.vcf >& tmp.log pigz -p 8 -f fluA.$asmAcc.$segRef.$today.nwk fluA.$asmAcc.$segRef.$today.vcf # Link to public trees archive directory (no assembly/segRef hierarchy, just by date) read y m d < <(echo $today | sed -re 's/-/ /g') archive=$archiveRoot/$y/$m/$d mkdir -p $archive ln -f $(pwd)/fluA.$asmAcc.$segRef.$today.{nwk,vcf,metadata.tsv}.gz $archive/ if [ -s fluA.$asmAcc.$segRef.$today.taxonium.jsonl.gz ]; then ln -f $(pwd)/fluA.$asmAcc.$segRef.$today.taxonium.jsonl.gz $archive/ fi gzip -c fluA.$asmAcc.$segRef.$today.pb > $archive/fluA.$asmAcc.$segRef.$today.pb.gz ln -f $(pwd)/hgPhyloPlace.description.$asmAcc.$segRef.txt \ $archive/fluA.$asmAcc.$segRef.$today.version.txt # Update 'latest' in $archiveRoot for f in $archive/fluA.$asmAcc.$segRef.$today.*; do latestF=$(echo $(basename $f) | sed -re 's/'$today'/latest/') ln -f $f $archiveRoot/$latestF done # Update hgdownload-test link for archive (adding assembly/segRef hierarchy) - mkdir -p $downloadsRoot/$asmDir/UShER_$segRef/$y/$m/$d - ln -sf $archive/*.$segRef.* $downloadsRoot/$asmDir/UShER_$segRef/$y/$m/$d/ - ln -sf $archiveRoot/*.$segRef.latest.* $downloadsRoot/$asmDir/UShER_$segRef/ + ssh hgwdev ln -sf $archiveRoot/*.$segRef.latest.* $downloadsRoot/$asmDir/UShER_$segRef/ # rsync to hgdownload hubs dir - for h in hgdownload1 hgdownload3; do - if rsync -a -L --delete $downloadsRoot/$asmDir/UShER_$segRef \ + for h in hgdownload1 hgdownload2 hgdownload3; do + if ssh hgwdev rsync -a -L --delete $downloadsRoot/$asmDir/UShER_$segRef \ qateam@$h:/mirrordata/hubs/$asmDir/; then true else echo "" echo "*** rsync to $h failed -- disk full ? ***" echo "" fi done done done -$fluAScriptDir/runGenoFlu.sh >& genoflu.log +$fluAScriptDir/runGenoFlu.sh $today >& genoflu.log echo "Built the trees, cleaning up." rm -f mutation-paths.txt *.pre*.pb final-tree.nh tmp.log nice gzip -f *.log *.tsv move_log* *.stderr echo "Building H5N1 outbreak trees" $fluAScriptDir/buildConcatTree.sh $today $fluAScriptDir/buildConcatTreeD1.1.sh $today echo "All done!"