e8b649b163d747c026dde4606542fcc024a6943a hiram Sun Sep 27 13:14:54 2026 -0700 add the mash sketch operation to the GenArk build refs #34360 diff --git src/hg/utils/automation/doAssemblyHub.pl src/hg/utils/automation/doAssemblyHub.pl index 883daa89349..6323b582f0f 100755 --- src/hg/utils/automation/doAssemblyHub.pl +++ src/hg/utils/automation/doAssemblyHub.pl @@ -7,54 +7,56 @@ # 1. Global-replace doTemplate.pl with your actual script name. # 2. Search for template and replace each instance with something appropriate. # Add steps and subroutines as needed. Other do*.pl or make*.pl may have # useful example code -- this is just a skeleton. use Getopt::Long; use File::Temp qw(tempfile); use File::stat; use warnings; use strict; use FindBin qw($Bin); use lib "$Bin"; use HgAutomate; use HgRemoteScript; use HgStepManager; +use AsmHub qw(accessionFromPath); # Option variable names, both common and peculiar to this script: use vars @HgAutomate::commonOptionVars; use vars @HgStepManager::optionVars; use vars qw/ $opt_buildDir $opt_sourceDir $opt_species $opt_rmskSpecies $opt_runRepeatModeler $opt_ncbiRmsk $opt_noRmsk $opt_augustusSpecies $opt_noAugustus $opt_xenoRefSeq $opt_noXenoRefSeq $opt_ucscNames $opt_dbName /; # Specify the steps supported with -continue / -stop: my $stepper = new HgStepManager( [ { name => 'download', func => \&doDownload }, { name => 'sequence', func => \&doSequence }, + { name => 'mashSketch', func => \&doMashSketch }, { name => 'assemblyGap', func => \&doAssemblyGap }, { name => 'chromAlias', func => \&doChromAlias }, { name => 'gatewayPage', func => \&doGatewayPage }, { name => 'cytoBand', func => \&doCytoBand }, { name => 'gc5Base', func => \&doGc5Base }, { name => 'repeatModeler', func => \&doRepeatModeler }, { name => 'repeatMasker', func => \&doRepeatMasker }, { name => 'simpleRepeat', func => \&doSimpleRepeat }, { name => 'allGaps', func => \&doAllGaps }, { name => 'idKeys', func => \&doIdKeys }, { name => 'windowMasker', func => \&doWindowMasker }, { name => 'addMask', func => \&doAddMask }, { name => 'gapOverlap', func => \&doGapOverlap }, { name => 'tandemDups', func => \&doTandemDups }, { name => 'cpgIslands', func => \&doCpgIslands }, @@ -134,30 +136,33 @@ expanded directory of mrnas/ and xenoRefMrna.sizes, default $xenoRefSeq _EOF_ ; print STDERR &HgAutomate::getCommonOptionHelp('dbHost' => $dbHost, 'workhorse' => $workhorse, 'fileServer' => $fileServer, 'bigClusterHub' => $bigClusterHub, 'smallClusterHub' => $smallClusterHub); print STDERR " Automates build of assembly hub. Steps: download: sets up sym link working hierarchy from already mirrored files from NCBI in: $sourceDir/GC[AF]/123/456/789/asmId sequence: establish AGP and 2bit file from NCBI directory + mashSketch: mash sketch the unmasked 2bit sequence into buildDir/mashSketch/, + the permanent cache AssemblyDivergence.pm/mashDistance.pl use to + pick an alignment pipeline/preset for a pair of assemblies assemblyGap: create assembly and gap bigBed files and indexes for assembly track names chromAlias: construct asmId.chromAlias.txt for alias name recognition gatewayPage: create html/asmId.description.html contents cytoBand: create cytoBand track and navigation ideogram gc5Base: create bigWig file for gc5Base track (and gcOnFly 2026-08-25) repeatModeler: optionally, run RepeatModeler to construct custom library for repeatMasker run. Use: -runRepeatModeler to perform this procedure, warning: can take a considerable amount of time (12 to 48 hours or more), and consumes an entire ku cluster node repeatMasker: run repeat masker cluster run and create bigBed files for the composite track categories of repeats simpleRepeat: run trf cluster run and create bigBed file for simple repeats allGaps: calculate all actual real gaps due to N's in sequence, can be @@ -987,30 +992,69 @@ exit 255 fi _EOF_ ); } $bossScript->add(<<_EOF_ twoBitToFa ../\$dbName.2bit stdout | faCount stdin | gzip -c > \$asmId.faCount.txt.gz touch -r ../\$dbName.2bit \$asmId.faCount.txt.gz zgrep -P "^total\t" \$asmId.faCount.txt.gz > \$asmId.faCount.signature.txt touch -r ../\$dbName.2bit \$asmId.faCount.signature.txt _EOF_ ); $bossScript->execute(); } # doSequence +######################################################################### +# * step: mashSketch [workhorse] +sub doMashSketch { + my $runDir = "$buildDir/mashSketch"; + &HgAutomate::mustMkdir($runDir); + + # AssemblyDivergence.pm/AsmHub.pm's mashSketchDir() expects the cache + # file named after the bare accession (e.g. GCA_939628115.1.msh), not + # the full asmId (e.g. GCA_939628115.1_Tfree1.0) -- accessionFromPath() + # only needs the accession prefix to match, so it works on the bare + # asmId string here the same way it works on a path's basename. + my $accession = &accessionFromPath($asmId); + if (! $accession) { + &HgAutomate::verbose(1, + "# mashSketch: '$asmId' is not a GCA/GCF accession, skipping\n"); + return; + } + + my $whatItDoes = +"mash sketch the unmasked 2bit sequence for use by AssemblyDivergence.pm/ +mashDistance.pl (mash distance -> alignment pipeline/preset triage) so it +never has to be built again for this assembly."; + my $bossScript = newBash HgRemoteScript("$runDir/doMashSketch.bash", + $workhorse, $runDir, $whatItDoes); + + $bossScript->add(<<_EOF_ +export asmId="$defaultName" +export accession="$accession" + +if [ ../\$asmId.unmasked.2bit -nt \$accession.msh ]; then + twoBitToFa ../\$asmId.unmasked.2bit stdout \\ + | mash sketch -k 21 -s 10000 -I \$accession -o \$accession - 2> /dev/null + touch -r ../\$asmId.unmasked.2bit \$accession.msh +fi +_EOF_ + ); + $bossScript->execute(); +} # doMashSketch + ######################################################################### # * step: assemblyGap [workhorse] sub doAssemblyGap { my $runDir = "$buildDir/trackData/assemblyGap"; &HgAutomate::mustMkdir($runDir); my $whatItDoes = "construct assembly and gap tracks from AGP file"; my $bossScript = newBash HgRemoteScript("$runDir/doAssemblyGap.bash", $workhorse, $runDir, $whatItDoes); $bossScript->add(<<_EOF_ export asmId="$defaultName" if [ ../../\$asmId.agp.gz -nt \$asmId.assembly.bb ]; then