1c41dff03108b4818777b9c9914e81c72555163f hiram Tue Sep 29 12:28:41 2026 -0700 miniMap2 now with swap option refs #34360 diff --git src/hg/utils/automation/doMiniMap2.pl src/hg/utils/automation/doMiniMap2.pl index 0ed60ca8612..8b14cec865f 100755 --- src/hg/utils/automation/doMiniMap2.pl +++ src/hg/utils/automation/doMiniMap2.pl @@ -34,30 +34,32 @@ use vars @HgAutomate::commonOptionVars; use vars @HgStepManager::optionVars; use vars qw/ $opt_buildDir $opt_target2Bit $opt_targetSizes $opt_query2Bit $opt_querySizes $opt_minimapPreset $opt_minimapCpu $opt_minimapSecondaryN $opt_minimapSecondaryRatio $opt_chainRam $opt_chainCpu $opt_regenerateMash + $opt_swap + $opt_swapDir /; # Specify the steps supported with -continue / -stop: my $stepper = new HgStepManager( [ { name => 'align', func => \&doAlign }, { name => 'chain', func => \&doChain }, { name => 'net', func => \&doNet }, { name => 'load', func => \&doLoad }, { name => 'cleanup', func => \&doCleanup }, ] ); # Option defaults: my $dbHost = 'hgwdev'; my $ramG = '24g'; # minimap2 index + 8 threads on a whole genome needs headroom @@ -123,49 +125,63 @@ axtChain/chainNet get the same kind of overlapping candidate alignments a lastz run would give them, and can resolve ambiguous/duplicated regions using genome-scale synteny context instead of relying on minimap2's own local primary/secondary call alone. -minimapSecondaryRatio F minimap2 -p: minimum secondary-to-primary score ratio to keep a secondary at all, default: $minimapSecondaryRatio -chainRam Ng Cluster ram size for chain step, default: -chainRam=$chainRam -chainCpu N Cluster CPUs number for chain step, default: -chainCpu=$chainCpu -regenerateMash Force the mash divergence check to re-sketch and overwrite even if a cached .msh (GenArk mashSketch/ cache or \$buildDir/mashDistance.txt) already exists. Ignored if -minimapPreset is also given, since that skips the mash run entirely. + -swap fromDb toDb are the same as the already-completed + primary run's; instead of aligning again, chainSwap + the primary run's axtChain/all.chain.gz (requires it + to already exist) and build the toDb-to-fromDb + liftOver/quickLift chain from that, in + -swapDir (default: toDb's own build root)/miniMap2.\$fromDb.swap + -swapDir dir Use dir instead of the default -swap work directory. _EOF_ ; print STDERR &HgAutomate::getCommonOptionHelp('dbHost' => $dbHost, 'workhorse' => '', 'fileServer' => '', 'ram' => $ramG, 'minimapCpu' => $minimapCpu, 'bigClusterHub' => ''); print STDERR " Automates a same-species/same-haplotype liftOver (minimap2/chain/net) pipeline, patterned after doSameSpeciesLiftOver.pl but using minimap2 in place of blat -fastMap: align: Aligns the assemblies using minimap2 -cx \$minimapPreset on a big cluster, one job per target sequence against the full query genome, then converts PAF+cigar to PSL with pafToPsl. chain: Chains the alignments on a big cluster. net: Nets the alignments, uses netChainSubset to extract liftOver chains. load: Installs liftOver chain files, calls hgAddLiftOverChain on $dbHost. cleanup: Removes or compresses intermediate files. +With -swap (not a normal stepper step -- see -swap above), chainSwap the +primary run's all.chain.gz and reuse the net/load logic to build the +reverse-direction liftOver/quickLift chain, without re-aligning. +Each completed run also leaves (or refreshes) a dateless symlink +miniMap2.\$ToDb -> the dated build directory (miniMap2.\$fromDb.swap for a +-swap run), so later runs/tools can find \"the\" build for a pair without +knowing its date. All operations are performed in the build directory, which defaults to $HgAutomate::clusterData/\$fromDb/$HgAutomate::trackBuild/miniMap2\$ToDb.\$date for a plain UCSC database \$fromDb, or to .../asmHubs/{genbankBuild,refseqBuild}/GCx/ddd/ddd/ddd/<asmId>/trackData/miniMap2\$ToDb.\$date when \$fromDb is a GenArk accession (matching how pairwise lastz builds are organized under GenArk) -- unless -buildDir is given. "; # Detailed help (-help): print STDERR " Assumptions: 1. /scratch/data/\$db/\$db.2bit contains sequence for database/assembly \$db. 2. $HgAutomate::clusterData/\$db/chrom.sizes contains all sequence names and sizes from \$db.2bit. 3. The \$db.2bit files have already been distributed to cluster-scratch (/scratch/data/<db>/). @@ -176,47 +192,49 @@ individual -- this is not a general any-vs-any pipeline. For that, use doBlastzChainNet.pl. " if ($detailed); print "\n"; exit $status; } # Globals: # Command line args: tDb=fromDb, qDb=toDb my ($tDb, $qDb); # Other: my ($buildDir); my ($tSeq, $tSizes, $qSeq, $qSizes, $QDb, $fileServer); -my ($liftOverChainDir, $liftOverChainFile, $liftOverChainPath); +my ($liftOverChainDir, $liftOverChainFile, $liftOverChainPath, $swapDir); sub checkOptions { # Make sure command line options are valid/supported. my $ok = GetOptions(@HgStepManager::optionSpec, 'buildDir=s', 'target2Bit=s', 'targetSizes=s', 'query2Bit=s', 'querySizes=s', 'minimapPreset=s', 'minimapCpu=i', 'minimapSecondaryN=i', 'minimapSecondaryRatio=f', 'chainRam=s', 'chainCpu=i', 'regenerateMash', + 'swap', + 'swapDir=s', @HgAutomate::commonOptionSpec, ); &usage(1) if (!$ok); &usage(0, 1) if ($opt_help); &HgAutomate::processCommonOptions(); my $err = $stepper->processOptions(); usage(1) if ($err); $dbHost = $opt_dbHost if ($opt_dbHost); if ($opt_minimapPreset) { $minimapPreset = $opt_minimapPreset; if ($minimapPreset !~ /^asm(5|10|20)$/) { die "-minimapPreset must be one of asm5, asm10, asm20 (got '$minimapPreset')\n"; } } # else: leave $minimapPreset undef -- estimateDivergence() will set it @@ -493,95 +511,104 @@ my $bossScript = newBash HgRemoteScript("$runDir/doChain.bash", $paraHub, $runDir, $whatItDoes); my $paraRun = &HgAutomate::paraRun($chainRam, $chainCpu); my $gensub2 = &HgAutomate::gensub2(); $bossScript->add(<<_EOF_ mkdir chainRaw $gensub2 pslParts.lst single gsub jobList $paraRun _EOF_ ); $bossScript->execute(); } # doChain +######################################################################### +# netFromChainBash($chainPath, $tSz, $qSz, $outDir, $overFile, $tName, $qName) +# -> bash source (string) that nets a single, already merge-sorted chain +# file into a liftOver chain + quickLift chain pair. $chainPath can be +# a literal path or a not-yet-Perl-interpolated bash expression (the +# caller controls that by how it quotes the argument) -- this doesn't +# care, it just emits it verbatim. Shared by doNet() (on the freshly +# merged chainRaw/*.chain) and doSwap() (on a chainSwap'd chain) so the +# "net one chain" logic isn't duplicated between them. +sub netFromChainBash { + my ($chainPath, $tSz, $qSz, $outDir, $overFile, $tName, $qName) = @_; + return <<_EOF_ +chainNet $chainPath \\ + $tSz $qSz \\ + noClass.net /dev/null +netChainSubset noClass.net $chainPath stdout \\ +| chainStitchId stdin stdout | gzip -c > $outDir/$overFile + +# make quickLift chain: +chainSwap $outDir/$overFile stdout \\ + | chainToBigChain stdin $outDir/$tName.$qName.quick.chain.txt \\ + $outDir/$tName.$qName.quick.link.txt +_EOF_ + ; +} # netFromChainBash + ######################################################################### # * step: net [workhorse] sub doNet { my $runDir = "$buildDir/axtChain"; my @outs = ("$runDir/$tDb.$qDb.all.chain.gz", "$runDir/noClass.net"); &HgAutomate::checkCleanSlate('net', 'load', @outs); &HgAutomate::checkExistsUnlessDebug('chain', 'net', "$runDir/run/chainRaw/"); my $whatItDoes = "It nets the chained minimap2 alignments and runs netChainSubset to produce liftOver chains."; my $mach = &HgAutomate::chooseWorkhorse(); my $bossScript = newBash HgRemoteScript("$runDir/doNet.bash", $mach, $runDir, $whatItDoes); - my $chromBased = (`wc -l < $tSizes` <= $HgAutomate::splitThreshold); - my $lump = $chromBased ? "" : "-lump=100"; $bossScript->add(<<_EOF_ unset TMPDIR if [ -d "/data/tmp" ]; then export TMPDIR="/data/tmp" elif [ -d "/scratch/tmp" ]; then export TMPDIR="/scratch/tmp" else tmpSz=`df --output=avail -k /tmp | tail -1` shmSz=`df --output=avail -k /dev/shm | tail -1` if [ "\$shmSz" -gt "\$tmpSz" ]; then mkdir -p /dev/shm/tmp chmod 777 /dev/shm/tmp export TMPDIR="/dev/shm/tmp" else export TMPDIR="/tmp" fi fi -# Use local scratch disk... this can be quite I/O intensive: tmpDir=`mktemp -d -p \$TMPDIR doMm2.net.XXXXXX` -# Merge up the hierarchy and assign unique chain IDs: -chainMergeSort run/chainRaw/*.chain \\ -| chainSplit $lump \$tmpDir/chainSplit stdin -endsInLf \$tmpDir/chainSplit/*.chain - -mkdir \$tmpDir/netSplit \$tmpDir/overSplit -for f in \$tmpDir/chainSplit/*.chain; do - split=\$(basename "\$f" .chain) - chainNet \$f \\ - $tSizes $qSizes \\ - \$tmpDir/netSplit/\$split.net /dev/null - netChainSubset \$tmpDir/netSplit/\$split.net \$f stdout \\ - | chainStitchId stdin \$tmpDir/overSplit/\$split.chain -done -endsInLf \$tmpDir/netSplit/*.net -endsInLf \$tmpDir/overSplit/*.chain - -cat \$tmpDir/chainSplit/*.chain | gzip -c > $tDb.$qDb.all.chain.gz -# noClass.net is scaffolding for netChainSubset above, not a kept product -# (matches doBlastzChainNet.pl's own treatment -- see cleanup step): -cat \$tmpDir/netSplit/*.net > noClass.net - -cat \$tmpDir/overSplit/*.chain | gzip -c > $runDir/$liftOverChainFile -# make quickLift chain: -chainSwap $runDir/$liftOverChainFile stdout \\ - | chainToBigChain stdin $runDir/$tDb.$qDb.quick.chain.txt \\ - $runDir/$tDb.$qDb.quick.link.txt +# Merge up the hierarchy and assign unique chain IDs -- one chain file for +# the whole genome. chainNet doesn't need this pre-split by chromosome +# the way doSameSpeciesLiftOver.pl split it; that was only ever a scale +# optimization doBlastzChainNet.pl itself only applies conditionally (its +# own \$splitRef), and this same-species/haplotype-scale pipeline doesn't +# need it either. +chainMergeSort run/chainRaw/*.chain > \$tmpDir/all.chain +gzip -c \$tmpDir/all.chain > $tDb.$qDb.all.chain.gz +_EOF_ + ); + $bossScript->add(&netFromChainBash('$tmpDir/all.chain', $tSizes, $qSizes, + $runDir, $liftOverChainFile, $tDb, $qDb)); + $bossScript->add(<<_EOF_ rm -rf \$tmpDir/ _EOF_ ); $bossScript->execute(); } # doNet ######################################################################### # * step: load [dbHost] sub doLoad { my $runDir = "$buildDir/axtChain"; &HgAutomate::checkExistsUnlessDebug('net', 'load', "$runDir/$liftOverChainFile"); my $whatItDoes = @@ -720,55 +747,164 @@ bedToBigBed -type=bed4+1 -as=bigLink.as -tab $tDb.$qDb.quick.link.txt $qSizes $tDb.$qDb.quickLink.bb qBases=`ave -col=2 $qSizes | grep "^total" | awk '{printf "%d", \$2}'` qCovered=`bigBedInfo $tDb.$qDb.quickLink.bb | grep "basesCovered" | cut -d' ' -f2 | tr -d ','` qPerCent=`echo \$qCovered \$qBases | awk '{printf "%.3f", 100.0*\$1/\$2}'` printf "%d bases of %d (%s%%) in intersection\\n" "\$qCovered" "\$qBases" "\$qPerCent" > $buildDir/fb.$tDb.quick${QDb}Link.txt rm -f bigChain.as bigLink.as $tDb.$qDb.quick.chain.txt $tDb.$qDb.quick.link.txt _EOF_ ); $bossScript->execute(); } # doLoad +######################################################################### +# swapGlobals: swap $tDb/$qDb (and everything derived from them) so +# doLoad() can be reused unchanged for the -swap direction, mirroring +# doBlastzChainNet.pl's own swapGlobals(). $buildDir becomes $swapDir, +# which the caller (doSwap()) must have already set. +sub swapGlobals { + ($tDb, $qDb) = ($qDb, $tDb); + $QDb = ucfirst($qDb); + ($tSeq, $qSeq) = ($qSeq, $tSeq); + ($tSizes, $qSizes) = ($qSizes, $tSizes); + $buildDir = $swapDir; + $liftOverChainDir = "$HgAutomate::clusterData/$tDb/$HgAutomate::trackBuild/liftOver"; + $liftOverChainFile = "${tDb}To${QDb}.over.chain.gz"; + $liftOverChainPath = "$liftOverChainDir/$liftOverChainFile"; +} # swapGlobals + +######################################################################### +# * -swap (not a normal stepper step): chainSwap the primary run's +# all.chain.gz and reuse doNet()'s netFromChainBash() + doLoad() to +# produce the reverse-direction liftOver/quickLift chain, without +# re-aligning anything. +sub doSwap { + my $origTDb = $tDb; + my $origQDb = $qDb; + + my $primaryChain = "$buildDir/axtChain/$origTDb.$origQDb.all.chain.gz"; + if (! -e $primaryChain) { + die "-swap: can't find $primaryChain -- run the primary $origTDb " . + "$origQDb build (through at least the 'net' step) first.\n"; + } + + $swapDir = $opt_swapDir ? $opt_swapDir : + &asmRoot($origQDb) . "/miniMap2.$origTDb.swap"; + my $runDir = "$swapDir/axtChain"; + &HgAutomate::mustMkdir($runDir); + + my $whatItDoes = +"It chainSwaps the primary $origTDb/$origQDb run's all.chain.gz and nets it +to produce the $origQDb-to-$origTDb liftOver/quickLift chain, without +re-aligning."; + my $mach = &HgAutomate::chooseWorkhorse(); + my $bossScript = newBash HgRemoteScript("$runDir/doSwap.bash", $mach, + $runDir, $whatItDoes); + + my $swappedOver = "${origQDb}To" . ucfirst($origTDb) . ".over.chain.gz"; + $bossScript->add(<<_EOF_ +unset TMPDIR +if [ -d "/data/tmp" ]; then + export TMPDIR="/data/tmp" +elif [ -d "/scratch/tmp" ]; then + export TMPDIR="/scratch/tmp" +else + tmpSz=`df --output=avail -k /tmp | tail -1` + shmSz=`df --output=avail -k /dev/shm | tail -1` + if [ "\$shmSz" -gt "\$tmpSz" ]; then + mkdir -p /dev/shm/tmp + chmod 777 /dev/shm/tmp + export TMPDIR="/dev/shm/tmp" + else + export TMPDIR="/tmp" + fi +fi +tmpDir=`mktemp -d -p \$TMPDIR doMm2.swap.XXXXXX` + +# chainSwap flips target/query coordinates on each already uniquely-ID'd +# chain block; chainSort re-sorts by the new target's coordinate order, +# same as doBlastzChainNet.pl's own swapChains() -- no need to re-run +# chainMergeSort, the IDs from the primary run are still unique. +chainSwap $primaryChain stdout | chainSort stdin \$tmpDir/all.chain +gzip -c \$tmpDir/all.chain > $origQDb.$origTDb.all.chain.gz + +_EOF_ + ); + $bossScript->add(&netFromChainBash('$tmpDir/all.chain', $qSizes, $tSizes, + $runDir, $swappedOver, $origQDb, $origTDb)); + $bossScript->add(<<_EOF_ +rm -rf \$tmpDir/ +_EOF_ + ); + $bossScript->execute(); + + # Reuse doLoad() unchanged for the swapped direction -- swapGlobals() + # makes it think it's doing a normal, primary $origQDb-to-$origTDb run: + &swapGlobals(); + &doLoad(); + + # Stable, dateless pointer to this swap build, in $origQDb's own root + # (relative target, so it stays correct if this tree gets mirrored): + my $swapBase = basename($swapDir); + my $stableLink = &asmRoot($origQDb) . "/miniMap2.$origTDb"; + &HgAutomate::run("rm -f $stableLink"); + &HgAutomate::run("ln -s $swapBase $stableLink"); +} # doSwap + + ######################################################################### # * step: cleanup [fileServer] sub doCleanup { my $runDir = "$buildDir"; my $whatItDoes = "It cleans up or compresses intermediate files."; $fileServer = &HgAutomate::chooseFileServer($runDir); my $bossScript = newBash HgRemoteScript("$runDir/doCleanup.bash", $fileServer, $runDir, $whatItDoes); $bossScript->add(<<_EOF_ rm -f run.mm2/query.fa rm -rf run.mm2/paf/ rm -rf axtChain/run/chainRaw/ rm -f axtChain/noClass.net # mashSketch.{a,b}.msh only ever exist here if estimateDivergence() fell # back to a one-off sketch (tSeq/qSeq didn't resolve to a real GenArk # accession or UCSC db -- see AssemblyDivergence.pm's sketch()); for a # normal assembly the sketch lives in its own permanent, shared cache # under /hive/data/genomes/, never here, so this is a no-op in that case: rm -f mashSketch.a.msh mashSketch.b.msh _EOF_ ); $bossScript->execute(); } # doCleanup +# asmRoot($db) -> this assembly's own build-directory root: +# a GenArk accession's real build dir + "/trackData" (matching how +# pairwise lastz builds are organized under GenArk), or +# clusterData/$db/$HgAutomate::trackBuild for a plain UCSC database. +# Used for both $buildDir (rooted at $tDb) and $swapDir (rooted at +# $qDb), and for where each run's dateless "current build" symlink goes. +sub asmRoot { + my ($db) = @_; + my $accession = &accessionFromPath($db); + my $asmBuild = $accession ? &asmBuildDir($accession) : undef; + return $asmBuild ? "$asmBuild/trackData" + : "$HgAutomate::clusterData/$db/$HgAutomate::trackBuild"; +} # asmRoot + # resolveAssemblySeq($db, $opt2Bit, $optSizes) -> ($seq, $sizes) # $opt2Bit/$optSizes, if given (-target2Bit etc.), always win. # Otherwise, if $db looks like a GenArk accession (accessionFromPath()) # with a real build directory (AsmHub::asmBuildDir()), resolve straight # to that build tree's own <asmId>.2bit/.chrom.sizes -- the same # resolution AssemblyDivergence.pm's sketch() already does for mash, # just for the actual alignment inputs here. Otherwise falls back to # the traditional /scratch/data or clusterData/$db/$db.2bit UCSC-db # location, as before. Doesn't require the result to actually exist -- # getSeqAndSizes() checks that afterward. sub resolveAssemblySeq { my ($db, $opt2Bit, $optSizes) = @_; my $accession = &accessionFromPath($db); my $asmBuild = $accession ? &asmBuildDir($accession) : undef; my $asmId = $asmBuild ? basename($asmBuild) : undef; @@ -886,54 +1022,62 @@ &getSeqAndSizes(); $QDb = ucfirst($qDb); $liftOverChainDir = "$HgAutomate::clusterData/$tDb/$HgAutomate::trackBuild/liftOver"; $liftOverChainFile = "${tDb}To${QDb}.over.chain.gz"; $liftOverChainPath = "$liftOverChainDir/$liftOverChainFile"; $chainRam = $opt_chainRam ? $opt_chainRam : $chainRam; $chainCpu = $opt_chainCpu ? $opt_chainCpu : $chainCpu; $minimapCpu = $opt_minimapCpu ? $opt_minimapCpu : $minimapCpu; $minimapSecondaryN = $opt_minimapSecondaryN ? $opt_minimapSecondaryN : $minimapSecondaryN; $minimapSecondaryRatio = $opt_minimapSecondaryRatio ? $opt_minimapSecondaryRatio : $minimapSecondaryRatio; $ramG = $opt_ram ? $opt_ram : $ramG; my $date = `date +%Y-%m-%d`; chomp $date; -if ($opt_buildDir) { - $buildDir = $opt_buildDir; -} else { - # GenArk target: build under its own trackData/, same convention - # ottoRequestWatch.sh expects pairwise lastz builds to already be in - # (".../<asmId>/trackData/lastz<Q>.YYYY-MM-DD") -- not under - # clusterData/$tDb/bed/, which only exists for native UCSC databases. - my $tAccession = &accessionFromPath($tDb); - my $tAsmBuild = $tAccession ? &asmBuildDir($tAccession) : undef; - $buildDir = $tAsmBuild ? - "$tAsmBuild/trackData/miniMap2${QDb}.$date" : - "$HgAutomate::clusterData/$tDb/$HgAutomate::trackBuild/miniMap2${QDb}.$date"; +$buildDir = $opt_buildDir ? $opt_buildDir : + &asmRoot($tDb) . "/miniMap2${QDb}.$date"; + +if ($opt_swap) { + # -swap reuses the primary run's already-completed axtChain/all.chain.gz + # (doSwap() checks for it directly) -- never auto-create $buildDir here. + if (! -d $buildDir) { + die "-swap: $buildDir does not exist -- run the primary $tDb $qDb " . + "build first, or give -buildDir to point at it.\n"; + } + &doSwap(); + &HgAutomate::verbose(1, + "\n *** All done! -swap steps were performed in $swapDir\n\n"); + exit 0; } if (! -d $buildDir) { if ($stepper->stepPrecedes('align', $stepper->getStartStep())) { die "$buildDir does not exist; try running again with -buildDir.\n"; } &HgAutomate::mustMkdir($buildDir); } &estimateDivergence(); $stepper->execute(); my $stopStep = $stepper->getStopStep(); my $upThrough = ($stopStep eq 'cleanup') ? "" : " (through the '$stopStep' step)"; +# Stable, dateless pointer to this build, refreshed on every run +# (relative target, so it stays correct if this tree gets mirrored): +my $stableLink = &asmRoot($tDb) . "/miniMap2.$QDb"; +&HgAutomate::run("rm -f $stableLink"); +&HgAutomate::run("ln -s " . basename($buildDir) . " $stableLink"); + &HgAutomate::verbose(1, "\n *** All done!$upThrough\n"); &HgAutomate::verbose(1, " *** Steps were performed in $buildDir\n"); if ($stepper->stepPrecedes('net', $stopStep)) { &HgAutomate::verbose(1, " *** Test installation ($HgAutomate::gbdb, goldenPath, hgLiftover " . "operation) on $dbHost.\n"); } &HgAutomate::verbose(1, "\n");